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Abstract 


A detailed description of a parameterization for thermal infrared radiative transfer 
designed specifically for use in global climate models is presented. The parameteri- 
zation includes the effects of the main absorbers of terrestrial radiation: water vapor, 
carbon dioxide, and ozone. While being computationally efficient, the schemes compute 
very accurately the clear-sky fluxes and cooling rates from the earth’s surface to 0.01 
mb. Errors are < 0.4 K day -1 in the cooling rate and < 1% in the fluxes. Since no 
transmittances are pre-computed, the atmospheric layers and the vertical distribution 
of the absorbers may be freely specified. The scheme can also account for any vertical 
distribution of fractional cloudiness with arbitrary optical thickness. In addition, the 
numerics and the FORTRAN implementation have been carefully designed to conserve 
both memory and computer time. This code should thus be particularly attractive 
to those contemplating long-term climate simulations, wishing to model the middle 
atmosphere, or planning to use a large number of levels in the vertical. 
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1 Introduction 


Thermal infrared radiation (IR) plays a crucial role in determining the earth’s climate and 
its sensitivity. It is, therefore, important to have an accurate IR radiation parameterization 
in atmospheric general circulation models (GCMs) used for studying climate change. In 
GCMs, calculations of thermal infrared fluxes can easily take a third or more of the com- 
puting time. As the spatial resolution and vertical extent of the models increase and the 
treatment of physical processes improves, it becomes imperative to have a fast and accurate 
IR radiation parameterization. 

Detailed calculation of the IR fluxes involves three integrations: a spectral integration, a 
vertical integration, and a directional integration. The spectral integration is the most 
time consuming. This is due to the narrowness of molecular absorption lines, which makes 
the absorption coefficient vary rapidly with wavenumber and thus requires very high spec- 
tral resolution 10 6 intervals for the entire IR spectrum) to obtain accurate integrated 
fluxes. But even this obstacle could be easily surmounted if it were not for the vertical 
integration. If the atmosphere were vertically homogeneous in pressure and temperature, 
the absorption coefficient over any vertical extent would be a function only of wavenumber, 
and wavenumbers with the same absorption coefficient would be radiatively identical. The 
spectral integration could then be greatly simplified by grouping all those spectral region 
with the same absorption coefficient — the ^-distribution approach (see, for example, Arking 
and Grossman 1972). In the real atmosphere, however, the dependence of the absorption co- 
efficient on pressure and temperature varies with wavenumber, and no two spectral intervals 
can be treated identically. It is this difficulty that makes the broad-band parameterization 
of IR fluxes such a difficult problem. 

The effect of vertical variations of pressure and temperature on absorption is commonly 
taken into account by using one of the following approximations: 


• One-parameter scaling , in which the effect of the variations of temperature and pres- 
sure along a path are taken into account by reducing a nonhomogeneous layer to 
an equivalent homogeneous layer with a “scaled” absorber amount at fixed reference 
temperature and pressure (e.g., Chou and Arking 1980). 

• Two-parameter scaling , in which the effects of variations of temperature and pressure 
along a path are taken into account by approximating a nonhomogeneous layer as a 
homogeneous layer with the actual absorber amount, but at an effective temperature 
and an effective pressure that depend on vertical variations within the layer (e.g., 
Kratz and Cess, 1988; Chou and Kouvaris 1991; Schwarzkopf and Fels 1991; Morcrette 
et al 1986; Kratz et al. 1993; Rosenfield 1991). In some cases (e.g., Schwarzkopf and 
Fels 1991; Rogers and Walshaw 1966), effective optical properties, such as line width 
and strength, are used as parameters instead of effective temperature and pressure; but 
the objective is the same — to account for vertical inhomogeneities in the atmosphere. 
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• Modified k- distribution methods , in which, as in the traditional ^-distribution ap- 
proach, wavenumbers with similar absorption coefficient are grouped, but in such a 
way as to account for the effects of vertical variations of temperature and pressure 
(Chou et al. 1993; Fu and Liou 1992; Goody et al. 1989; Lacis and Oinas 1991; Wang 
and Shi 1988; Zhu 1992). 


At the Goddard Climate and Radiation Branch, we have developed various broad-band IR 
radiation parameterizations for the major water vapor, carbon dioxide, and ozone absorption 
bands based on all three approximations (Chou and Kouvaris 1991; Chou et al. 1993; Chou 
et al. 1994). These IR radiation parameterizations have been shown to be accurate and 
efficient in computing cooling rates in the troposphere and lower stratosphere, as well as in 
the middle atmosphere (up to the 0.01-mb level). 

This report describes in detail how these parameterizations have been integrated into a 
single IR radiation package for use in GCMs. This package has been implemented and 
tested in the Goddard GCMs and is currently being used for both climate modeling and 
data assimilation applications. In Section 2 we briefly review the transfer equations for fluxes 
and cooling rates in clear and cloudy atmospheres. Section 3 addresses the partitioning of 
the spectrum into broad bands and the representations of the spectrally integrated fluxes in 
these bands. The forms of the transmittances used and the scheme for overlapping gaseous 
absorptions are given in Section 4. The vertical discretization of the transfer equations is 
given in Section 5. Section 6 presents some results of flux and cooling rate calculations and 
compares the parameterization against line-by-line calculations. The results are summarized 
in Section 7. 


2 Infrared Transfer Equations 

2.1 Radiative Transfer Through a Clear Atmosphere 

Let us consider a thin atmospheric layer at pressure p f with a temperature T", which has a 
differential pressure thickness dp f and an optical thickness du u at wavenumber v (Figure 1). 
The radiance emitted by this layer is , where p is the cosine of the angle between 

the beam and the vertical, and B u (T f ) is the Planck function. Let us further consider a 
higher level at pressure p where we wish to compute the upwelling radiance, and let us 
assume that between p and p ! there is a mixture of gaseous absorbers with monochromatic 
optical thickness 

t 


2 


lir 


( 2 ) 


where ul t (p,p f ) is the optical thickness of the zth absorber given by 

Jp 9 

Here g* is the zth absorber’s specific mass, k Kt is its absorption coefficient, and g is the 
gravitational acceleration. The quantity is the differential atmospheric mass per unit 
horizontal area. 



P 


P’ 

P’ + dp’ 


Figure 1: Contribution from a differential layer at p ' to the upward radiance at p. 


The transmittance between p and p f in the direction /z is then e and the contri- 

bution of the layer between p ' and p l + dp ' to the upwelling radiance at p is given by 

= ff„(r')— e“ u?(p ’ p ' )/M . (3) 

P 

The total contribution of the layer to the upward flux at p can be obtained by angular 
integration of (3) over the hemisphere, 

dF]{p ) = 27 r / dll(fj.)fidfi, (4a) 

Jo 

= 2'kB„{T') f 1 (4b) 

Jo P 

= -B v (T , )dT v (p,p '), (4c) 

where B U (T ) = and r^(p, p f ) is the flux transmittance for isotropic radiation given 

by i 

r v(p,P , ) = 2 [ e~ u ^ Plp, ^^pdfx. (5) 

Jo 

Equations analogous to (3)-(5) can be written for the downward radiance and flux at p due 
to differential layers above p. 
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2.2 The Effects of Clouds 


To model the effects of clouds we follow the general treatment described in Harshvardhan et 
al. (1987). Clouds are allowed to occupy any combination of model layers, and the fractional 
coverage and cloud optical thickness are assumed to be specified for each layer. 

Eq (5) applies to the case of a clear atmosphere, but it can be easily extended to the case of 
a cloudy and horizontally homogeneous (i.e., overcast) atmosphere by replacing the gaseous 
optical thickness with the total optical thickness: 

u v {p, p') = u 9 v (p, p') + ul(p, p'), (6) 

where 

3 

is the total cloud optical thickness between p and p\ and the u^ / j(p } p f ) is the optical 
thicknesses of the jth cloud layer. 

If, in addition to clear and overcast layers, a single cloud layer with fractional cover A c and 
optical thickness u c u is introduced between p and p’ } the mean radiance transmittance in 
the direction /z becomes 

(l _ c -u*,(p,p')//u _|_ i 4 C e -[u l/ (p,p ; )+u=]/M^ 

and the mean flux transmittance analogous to (5) is 

T v(P> P ') = 2 /' [(i - A c ) + ^ e -K(p.p')+^l/p] ^ 

= [l - A c (l - e _ “- /M )j T u (p,p'), (7) 

where JI is the effective mean value of p that converts the radiance transmittance to flux 
transmittance. Clearly, the diffusivity factor = depends on both u v and u c v , but it is com- 
monly taken to be constant and assigned the value 1.66. 

Letting N u = A c { 1 - (7) becomes 

K(p,p') = (1 - N v )t v {p,p'). (8) 

It is important to note that in this case the cloud fraction, A c , and the cloud optical 
thickness, u c v , enter the equations only through the quantity N v . Since the factor ( l-e _U| ' /**) 
is the flux emissivity of the cloud layer, N„ can be regarded as the cloud fraction of an 
equivalent black cloud; alternatively, (1 - JV v ) can be regarded as the effective transmissivity 
of an equivalent overcast cloud. 
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When there is more than one cloud layer with fractional cover between p and p\ the situation 
is considerably more complicated, since we need to describe how the clouds are overlapped. 
In general, we can write 

t Z(p*p') = C l/ (p,p f )r u (p i p , ) > (9) 

where C u {p^ p*) depends on the values of A c and u c v for the various cloud layers. 

For the special case of any combination of overcast and randomly overlapped fractional cloud 
layers, we have 

= IK 1 - **>). (io) 

J 

where the subscript j denotes the cloud layers. 

Since Y\j(l — Nvj) is the fraction of the horizontal area that would be cloud-free if all cloud 
layers (including overcast layers) were assigned their equivalent black-cloud fraction and 
then randomly overlapped, Harshvardhan et al. (1987) referred to C u {p^p f ) as the probability 
of a clear line of sight between p and p f . The clear-line-of-sight interpretation of C u {p,p ! ) is 
also valid for any overlapping of black clouds, making C u {p, p ! ) easy to obtained. Obtaining 
CvijPiP 1 ) for non-random arrangements of gray clouds, however, is not as straightforward. 
In general, C^(p,p') does not depends on N u alone, but on more complicated combinations 
of the A c - and the u c u -. 

As an example, consider another arrangement discussed by Harshvardhan et al. (1987), 
that of maximally overlapped clouds. If the largest cloud between p and p f is black, then 
^(PjPO is simply equal to one minus its fraction. If all clouds are allowed to be gray, 
however, C u {p,p f ) will depend on all cloud fractions and optical thicknesses. For a case 
with J maximally overlapped clouds between p and p', C u (p,p f ) may be easily computed 
by putting the clouds in order of increasing cloud cover and then evaluating the following 
recursion: 

JVj 0) = 0, (11a) 

A[W = Nl/j +ArU-i) e -<,i/* } j = i,...,j. (11b) 

One can easily verify that 

C„(p,p')= ( 12 ) 

Note that whenever a black cloud —> oo) occurs in (11) all smaller clouds are eliminated 
from the recursion. 

In the current code we have only implemented the random overlap strategy; but the com- 
putation of C(p,p f ) is well isolated so that it can be easily replaced if maximal overlapping 
or some other arrangement is desired. 
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2.3 The Flux Calculation 


Using the cloudy transmittance r* from (9) in (4c), the contribution of the layer between 
p f and p l + dp ' to the upward flux at p becomes: 

dFj(p) = -B v {T')dr :{ P , P '). (13) 

Similarly, it can be easily shown that the contribution of the earth’s surface to the upward 
flux at level p is given by 

6F}(p) = B v (T a )T*(p,p,), (14) 

where the subscript s denotes surface quantities. The total upward flux at p is then obtained 
by integrating over all wavenumbers and all differential layers below p: 

FHp) = jdu [B v (T t )r:(p,p t ) - jT' B u (T')^^-dp'j . (15a) 

Similarly, assuming no downward flux at p — 0, the total downward flux at p is 

F > {p) = J d ^J’ B „(T') e -?^dp}. (15b) 

The cooling rate is simply proportional to the flux divergence: 

= (16) 

where c p is the heat capacity of air at constant pressure. 


3 Broad-Band Calculations 

3.1 The 8 Bands 

For computing IR fluxes due to water vapor, carbon dioxide, and ozone, we divide the 
spectrum into eight bands. Table 1 shows the spectral ranges of these bands, together with 
the absorbers involved in each band and the parameterization methods used to compute the 
transmittance in each band. The transmittance parameterizations are discussed in Section 
4. 

Water vapor line absorption covers the entire IR spectrum, while water vapor continuum 
absorption is included in the 540-1380 cm -1 spectral region (Bands 3 through 6). The 
absorption due to C0 2 is included in the 540-800 cm' 1 region (Band 3), and the absorption 
due to O3 is included in the 980-1100 cm -1 region (Band 5). The division of Band 3 into 
three sub-bands is discussed in Section 4.5. 
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Table 1: Spectral bands, absorbers, and parameterization methods. 



Spectral 


Option for 

Band 

Range 

Absorber 

Transmittance 


(cm" 1 ) 


Parameterization 




"LOW” 

“HIGH” 

i 

0-340 

H 2 O line 

K 

T 

2 

340-540 

H 2 O line 

K 

T 

3a 

540-620 

^ H 2 0 line 

K 

K 

3b 

620-720 

> H 2 O continuum 

S 

S 

3c 

720-800 

J co 2 

K 

T 

4 

800-980 

H 2 0 line 

K 

K 



H 2 0 continuum 

S 

S 



H 2 O line 

K 

K 

5 

980-1100 

H 2 O continuum 

S 

S 



0 3 

T 

T 

6 

1100-1380 

H 2 O line 

K 

K 



H 2 O continuum 

S 

S 

7 

1380-1900 

H 2 0 line 

K 

T 

8 

1900-3000 

H 2 0 line 

K 

K 

K: ^-distribution method with linear pressure scaling. 


T: Table look-up with temperature and pressure scaling. 


S: One-parameter temperature scaling. 
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3.2 Planck- Weighted Band Integrals 


From (15a, b), the integrated fluxes for each of the bands can be written as: 
Jf0>) = B,(T.)r;(p,p.) - f B,(T') f 9r j ) P ' r * ) dp', 

where T' is the temperature at p', 

F t (p) = [ F v (p)du, 

Bi{T) = [ B v {T)dv y 

J Ai/; 

Ti(P.P) - flpi) 

fdT , (p,p’)\ _ / a^ F u (T') dT ^p p ^du 

V v /, ” 5 -( T ') 

are the Planck-weighted band integrals, and A is the width of the ith band. 


(17a) 

(17b) 

(18a) 

(18b) 

(18c) 

(18d) 


Within each spectral band, either the range of B u is sufficiently small or the shape of B v is 
sufficiently independent of temperature that we can make the following approximations: 


T*( P ,P') 


f dr*{p,p') 
\ dp' 


/ a , B u (T 0 )T*(p,p')dv 
B^To) 

Ia^ By(T D ) aT| 3 p/ P ^ dt/ _ frr*(p,p') 
5,(T 0 ) 5p' 


(19a) 

(19b) 


where T c is a typical value of the atmospheric temperatures. (In parameterizing r, below, 
we use T 0 = 250 K.) Figure 2 illustrates the accuracy of this approximation for water 
vapor line absorption. It compares the band-integrated flux transmittance of Bands 1, 2, 
and 3 weighted by £„(250K) with the same transmittances but weighted by fl„(200K) and 
5^(300K). The calculations were done at pressures of 10 mb and 500 mb and with the 
water vapor amount varying from 10~ 6 g cm -2 to 10 g cm -2 . Clearly, (19a) is an excellent 
approximation for a wide range of temperatures, pressures, and humidities. 


With the approximations (19a,b), the band-integrated fluxes reduce to: 


fHp) 
fRp ) 


^(r,K(p,p,) - J*‘ B i (T') dr: ]*& dj 


(20a) 

(20b) 
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0.0 0.2 0.4 0.6 0.8 1.0 

Flux Transmittance Weighted by B(250 K) 


Figure 2: Relation between flux transmittances weighted by B l/ (250K) to those weighted by 

B^SOOK) and B*,(200K). Results plotted are for various pressures and water vapor amounts. 


The parameterization problem is thus reduced to obtaining the T*(p,p 7 ) for each of the eight 
bands. 


Unlike the gaseous absorption coefficient, the cloud absorption coefficient varies slowly with 
wavenumber, and the cloud emissivity can be considered constant within a band. Using (9) 
and the approximation (19a), the flux transmittance becomes 


A B u (T 0 )C u {p,p')T u (p,p')dv 

Tl[P,P) ~ B X (T 0 ) 

« 5 = C t (p,p)r t (p,p). 


( 21 ) 


Derivations of tv, the Planck-weighted flux transmittance due to gaseous absorption, are 
given in Section 4. 
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3.3 Band-Integrated Planck Functions 

The spectrally integrated Planck fluxes were pre-computed for each band and then fitted 
by a 4th-degree polynomial in temperature: 

Bi(T) = c ti o + £ c l>n T n . (22) 

71=1 

When integrated over the eight bands, errors in this regression are negligible (< 0.1%) for 
160 K < T < 345 K. The coefficients c l>n are listed in Table 2. 


Table 2: Coefficients for computing the spectrally integrated Planck fluxes using the polynomial 
fits. The units of temperature are Kelvin. 


Band 

c,,o 

C»,l 

C.,2 

C », 3 

c i t 4 

i 

-2.6844E-1 

-8.8994E-2 

1.5676E-3 

-2.9349E-6 

-2.2233E-9 

2 

3.7315E+1 

-7.4758E-1 

4.6151E-3 

— 6.3260E-6 

3.5647E— 9 

3 

3.7187E+1 

-3.9085E-1 

— 6.1072E-4 

1.4534E-5 

— 1.6863E-8 

4 

— 4.1928E+1 

1.0027E+0 

-8.5789E-3 

2.9199E-5 

-2.5654E-8 

5 

-4.9163E+1 

9.8457E-1 

— 7.0968E-3 

2.0478E-5 

— 1.5514E— 8 

6 

-1.0345E+2 

1.8636E+0 

— 1.1753E— 2 

2.7864E— 5 

— 1.1998E-8 

7 

-6.9233E+0 

-1.5878E-1 

3.9160E-3 

— 2.4496E-5 

4.9301E-8 

8 

1.1483E+2 

-2.2376E+0 

1.6394E-2 

-5.3672E-5 

6.6456E-8 


4 Transmission Functions 


For accuracy and speed, the Planck- weighted flux transmittance for gaseous absorption: 


'r.(p.p') 


B u (T 0 )T v {p,f/)dv 
Bi{ T 0 ) 


(23) 


are computed using three different approaches, depending on the absorber and the spectral 
band: 
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• The k- distribution method with linear pressure scaling is applied to the water vapor 
bands. This method is also applied to the 15 fim CO2 band if accurate cooling rate 
calculations in the middle atmosphere are not required. 

• The transmittances due to CO2 and O3 absorption in Bands 3 and 5 are obtained 
from pre-computed transmittance tables based on two-parameter scaling . Because 
the ^-distribution method with linear pressure scaling underestimates the water vapor 
cooling rate in the middle atmosphere, the transmittances of the three strongest water 
vapor absorption bands (Bands 1, 2, and 7) are also obtained from pre-computed 
transmittance tables if accurate computations of the water vapor cooling in the middle 
atmosphere are required. 

• The transmittances due to water vapor continuum absorption in Bands 3 through 6 
are computed using a one-parameter scaling approach. 


The applications of these parameterizations to the different spectral bands and absorbers 
are summarized in Table 1. The code allows “HIGH” and “LOW” options to be specified 
depending on the desired accuracy in the middle atmosphere. With the “HIGH” option 
the more accurate, but more expensive, table look-up method is applied to C0 2 absorption 
and to water vapor line absorption in Bands 1, 2, and 7. With the “LOW” option these 
are all done with the ^-distribution method. The intermediate option of doing only C0 2 
by table look-up, which may be adequate for some stratospheric modeling when economy 
is important, is not explicitly provided but could be easily implemented. 


4.1 The ^-Distribution Method 

As shown in Chou et al. (1993), the cooling due to water vapor in the lower atmosphere 
(p > 20 mb) is primarily attributable to the spectral regions away from the center of 
absorption lines — where the absorption coefficient is approximately linear in pressure and 
its dependence on temperature varies smoothly with wavenumber. Under these conditions, 
the absorption coefficient at any temperature and pressure is simply proportional to its 
value at a reference pressure and temperature. This scaling of the absorption coefficient is 
also appropriate for C0 2 absorption in the troposphere and lower stratosphere. 

When this scaling is applicable, we adopt the following form for the absorption coefficient: 

k u (p,T) = k u (p r ,T T ) (^) m h(T,T r) , (24) 

where p T and T r are the reference pressure and temperature, m is an empirical constant, 
and h(T,T r ) is the temperature scaling factor, which satisfies h(T r ,T r ) = 1 . We use the 
following empirical relation for h(T,T T ): 

h(T,T r ) = 1 + a[T — T r ) + /3(T — T r ) 2 , (25) 
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where the coefficients a and 0 — constants for each absorber and spectral band— are obtained 
by minimizing errors in the flux transmittances. Details of their derivation are given in Chou 
et al. (1993). Values of p r , T r , m, a, and (3 used to scale the absorption coefficients for 
water vapor and CO2 in each of the bands are given in Table 3. Note that Band 3 is divided 
into three sub-bands; the reasons for this are explained in Section 4.5. With the scaling of 


Table 3: Scaling parameters used with the ^-distribution method. Bands are defined in Table 1. 
The reference temperature and pressure, T r and p r , are those used in (24), m is the power for the 
pressure scaling in the same equation, and a and 0 are the temperature scaling coefficients that 
appear in (25). 



h 2 o 

p T — 500 mb 
T r = 250 K 
m = 1 

C0 2 

T t = 250 K 


a (K- 1 ) 

P (K- 2 ) 

p T (mb) m 

a (K _1 ) 

P{ K- 2 ) 

Band 






1 

.0021 

— 1. Ole-5 




2 

.0140 

5.57e-5 




3a 



300 .50 

.0182 

1.07e-4 

3b 

\ .0167 

8.54e— 5 

30 .85 

.0042 

2.00e-5 

3c 

J 


300 .50 

.0182 

1.07e— 4 

4 

.0302 

2.96e— 4 




5 

.0307 

2.86e— 4 




6 

.0154 

7.53e-5 




7 

.0008 

— 3.52e-6 




8 

.0096 

1.64e— 5 





(24), the monochromatic flux transmittance given by (5) reduces to 

r„(p, p') = 2 C e - fc ^(P. P ') (26) 
Jo 

where q(p,p') is the scaled absorber amount between p and p' given by 

^y)=/V)(9V(Ar,)f, m 

and k T v = k v (p r ,T r ). In (27), we have assumed a single absorber. The overlapping of gaseous 
absorbers is discussed in Section 4.4. 
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It can be seen from (26) and (27) that by scaling the absorption coefficient the dependence 
of r u on wavenumber (through k T v ) is separated from its dependence on pressure and tem- 
perature (through q). Wavenumbers which have the same absorption coefficient at p r and 
T r will thus have the same transmittance at all other pressures and temperatures. Within 
spectral intervals narrow enough to consider the Planck function constant, these wavenum- 
bers are radiatively identical, and the integration over wavenumber can be replaced by an 
integration over the absorption coefficient k at the reference temperature and pressure: 


f [%)du — Av f ( 9)f r (k)dk . 

J Av J 0 


(28) 


Here f r (k)dk is the the fraction of the spectrum in band Au with absorption coefficients 
between k and k + dk , so that 

yoo 

f T {k)dk = 1. (29) 


f 


The function f T (k ) is called the ^-distribution function. We have used the superscript r as 
a reminder that we are using the ^-distribution at the reference temperature and pressure. 
This is the basic formalism of the ^-distribution method. 

Unfortunately, our 8 bands are too broad to assume a constant Planck function, and so 
(28) cannot be applied directly. Instead, we divide each band into a number of narrow 
sub-bands, for which the Planck function can be considered constant, and apply the k- 
distribution method to each sub-band. Letting 5vj be the width of each sub-band, the 
Planck-weighted transmission function given in (23) becomes 


, ffo. r v{q)dv 

r(,) = WT) ’ 


(30) 


where we have dropped the band index and (30) is understood to apply to each of the 
eight main bands. In (30), Bj(T 0 ) is the constant reference Planck flux in sub-band j , and 
B(T 0 ) — Yjj Bj(T 0 )Si/j is the band-integrated reference Planck flux. 


Using (28), the spectral integral is replaced with an integration over all values of k\ 

'M r A > 


-(?') = £ 


B(T 0 ) 


where 


T k{q) 


yoo 

fiVj / Tk(Q)fj{ k ) dk 
Jo J 

! f 1 

Jo 


(31) 


We use a different ^-distribution function, fj(k) } in each sub-band, but the same scaling of 
the absorber for all sub-bands within each band. 

As shown in Chou and Arking (1980), the flux transmittance can be approximated by 


r k {q) 


0 -kq/n 


(32) 
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where = is the diffusivity factor taken to be 1.66. Using this approximation, (31) reduces 


to 


r (?) = E 

3 




(33) 


The integral over k can then be replaced by a simple quadrature over N ^-intervals of width 
(6k) n located at k n , where n = 1, . . N. In doing this we use the same fc-intervals for all 
sub-bands within a band, but allow the intervals to vary from band to band. The band 
transmittance can then be expressed as an exponential sum: 


T(g) = X>T WM (A<7)n, (34) 

n=l 

where (A g) n is the Planck-weighted ^-distribution function for the nth fc-interval given by 

(A S )„ = E/J(W(«)n^^. (35) 

Note that for each fc-interval the (A^f) n are summations over all sub-bands and are indepen- 
dent of temperature, pressure, and absorber amount; they can therefore be pre-computed 
for each band. 


Values of (A g) n were first obtained from line-by-line calculations and then slightly adjusted 
using regression so that the rms difference between the transmittances computed from the 
line-by-line method and from (34) is minimized. The adjustment was necessary because 
of the diffusivity approximation (= = 1.66) and the small number of terms ( N < 6) used. 
Details are given in Chou et al. (1993). 

Calculations of the Planck- weighted flux transmission function using (34) can be done very 
efficiently, particularly if the k n are appropriately chosen. If the transmittances need to be 
computed between L + 1 levels separating L contiguous layers, (34) involves \ L(L - 1 )N 
exponentials. Note, however, that for multi-layer regions the exponentials in each term 
in (34) can be obtained by multiplying together the exponential factors for the individual 
layers; thus, only LN exponential calculations are required. To reduce this still further, we 
choose the k n so that they satisfy 

k n = rjk n _i , n = 2, N, (36) 

where rj is a positive integer. With this choice, only a single set of L exponential operations 
for k\ is needed. The other exponential terms can be obtained by raising the first to an 
integer power. For the values of tj used, this can be done with only 3 or 4 multiplication for 
each succeeding exponential. This greatly accelerates the transmittance calculation. 

In the parameterization, we apply the /^-distribution method as just described to compute 
the water vapor fine absorption in all bands but Band 3. The first absorption coefficient, 
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&i, the constant 77 , and the (A^) n are given in Table 4. As shown in Table 1, in Bands 
1, 2, and 7 water vapor line absorption can be computed using either this method or the 
transmittance tables described in the next subsection; in all other bands we use only the 
^-distribution method. In Band 3, the 15 fx m band, we use a somewhat more complicated 
procedure described in Section 4.5. 


Table 4: Parameters used for the transmission functions due to water vapor line absorption. 4* is 
the first flux absorption coefficient used with the ^-distribution method. The remaining coefficients 
are defined using k n = 77& n -i> where 77 is the band dependent constant given in the table. A g is 
the Planck-weighted ^-distribution function defined in (35). The use of the ^-distribution method 
in Band 3 is described in Section 4.5. Units of k are g -1 cm 2 ; other values are non-dimensional. 


Band 1 Band 2 Band 4 Band 5 Band 6 Band 7 Band 8 


M 

2.96e+l 

4.17e— 1 

5.25e— 4 

V 

6 

6 

6 

&g 1 

.2747 

.1521 

.4654 

&92 

.2717 

.3974 

.2991 

&g 3 

.2752 

.1778 

.1343 

a# 4 

.1177 

.1826 

.0646 

A55 

.0352 

.0374 

.0226 

A<?6 

.0255 

.0527 

.0140 


5.25e-4 

2.34e— 3 

1.32 

5.25e-4 

6 

8 

6 

16 

.5543 

.1846 

.0740 

.1437 

.2723 

.2732 

.1636 

.2197 

.1131 

.2353 

.4174 

.3185 

.0443 

.1613 

.1783 

.2351 

.0160 

.1146 

.1101 

.0647 


.0310 

.0566 

.0183 


4.2 Pre-Computed Transmittance Tables 

Although the ^-distribution method with linear pressure scaling is computationally very 
fast, it is not accurate in the middle atmosphere (0.01-20 mb), where the pressure ranges 
by three orders of magnitude and where the Doppler broadening of absorption lines is 
important. Radiative cooling in the middle atmosphere is primarily due to CO2 in the 15 
/zm band (Band 3), to O3 in the 9.6 fxm band (Band 5), and to a lesser extent, to water 
vapor near the centers of absorption bands (Bands 1, 2, and 7). 

As shown in Chou and Kouvaris (1991) and Chou et al. (1994), the transmittances in these 
bands can be simply derived from pre-computed transmittance tables. In this approach the 
tables are based on two-parameter scaling in which a non-homogeneous layer with pressure 
and temperature varying with height can be treated as if it were homogeneous with an 
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effective pressure and temperature given by 


( 37 ) 


Peff 


T eff 


J pdw 
Jdw ’ 
[Tdw 
Jdw ’ 


(38) 


where w is the absorber amount, and the integration is over the depth of a layer. It has 
been well recognized (e.g., Wu 1980; Chou and Kouvaris 1991) that the cooling rates depend 
primarily on contributions to the fluxes from nearby layers. Since pressure and temperature 
variations over nearby layers are generally small, the simple scaling approximations (37) and 
(38) produce accurate cooling rates. 


With two-parameter scaling, the flux transmittance becomes a function of the absorber 
amount and the effective temperature and pressure. These dependencies can be accurately 
pre-computed from the following equation using a line-by-line method, 


r(w,p c ({,T e {[) = 


_ jAv T v{' W >P'H,T'f{)B l ,{ T 0 ) d V 


Ta.. BjT n )dv 


(39) 


where Ay is the width of the entire spectral band. The band transmittance varies rapidly 
with pressure, but rather smoothly with temperature; thus, the size of the three-dimensional 
transmittance tables can be reduced to three two-dimensional tables using a quadratic fit 
in temperature, 

t(w,PcE, T e fr) = a(w, p e ff) + b(w,p c s)(T e ff - 250 K) + c(u;,p eff )(T eff - 250 K)\ (40) 

where a, b, and c are regression coefficients. This regression is valid for temperatures ranging 
from 170 K to 330 K; for this range it introduces only « 1% error m the absorptance. 
Because of the wide range of variation of w and p e ff, values of the coefficients are tabled at 
equal intervals in log 10 (w) and log 10 (p e ff)- Table 5 specifies the ranges used for the tables 
for C0 2 absorption in the 15 pm region, 0 3 absorption in the 9.6 pm region, and water 
vapor absorption in the three strong absorption bands. The actual tables can be obtained 
from the authors. 


This method is simple and accurate. As shown in Chou and Kouvaris (1991) and Chou et 
al. (1994), it can be used to calculate fluxes and cooling rate in both the middle and lower 
atmosphere, from 0.01 mb to the earth’s surface. However, this method can be significantly 
slower than the fc-distribution method on computers that cannot efficiently perform table 
look-ups; we have therefore provided the option of using the ^-distribution method for all 
bands, except for the computation of ozone transmittance in Band 5. 


4.3 One-Parameter Scaling For Water Vapor Continuum Absorption 


As shown in Table 1, water vapor continuum absorption is included in Bands 3 through 6. 
The water vapor continuum absorption coefficient, depends on the water vapor partial 
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Table 5: The first values, intervals, and sizes of the pre-computed tables of the regression coefficients 
a, 6 , and c used in (40) for transmittance calculations. Units of iy e ff are g cm -2 for water vapor and 
(cm-atm)sTP for CO 2 and O 3 ; p e ff is in mb. 



Band 1 

Band 2 

Band 3 

Band 5 

Band 7 

ABSORBER 

h 2 o 

h 2 o 

co 2 

O3 

h 2 o 

iogioC^eff)! 

-8 

-8 

-4 

-6 

-8 


-2 

-2 

-2 

-2 

-2 

Alogio(™ cff ) 

0.3 

0.3 

0.3 

0.3 

0.3 

Alog 10 (p eff ) 

0.2 

0.2 

0.2 

0.2 

0.2 

^-dimension 

31 

31 

24 

21 

31 

p-dimension 

26 

26 

26 

26 

26 


pressure, p e , and on temperature. It increases with increasing partial pressure and with 
decreasing temperature. As shown by Roberts et al. (1976), it can be approximated by 


K(Pe> T ) = Kfl — ex P 

Vo 


1800 


f- - —)] 

\T 296 /J ’ 


(41) 


where A£ 0 is the absorption coefficient when p e is p Q ( = 1013 mb), and the temperature is 
296 K. Values of k^ 0 are computed from the analytical representation given by Roberts et 
al. (1976), which is a fit to the laboratory data of D. E. Burch. 


With the scaling of (41) for the absorption coefficient, the flux transmittance reduces to 

t u {p, p') = 2 [ l e -^,o XP’P'Vvpdp, (42) 

Jo 

where q is the scaled water vapor amount for continuum absorption given by 

«(p,p') = — f V'. 

9 Jv Vo 


g J P Vo 

Here we have used the approximation 


where q is the specific humidity. 


q(p)p 

0.622’ 


(43) 


(44) 
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( 45 ) 


Taking the Planck-weighted average of (42) over a band we have 

_ Iau T v (q)B v (T 0 )dv 
T[q) ~ J&v B v (To)dv * 

where we note that the transmittance is a function of q only. The absorption coefficient 
varies slowly with wavenumber, except in Band 3 where it varies by a factor of 3.5; we there- 
fore divide Band 3 into three sub-bands, as described in Section 4.5. The Planck-weighted 
flux transmittances for the three sub-bands and for Bands 4, 5, and 6 were computed from 
(45). These transmittances were then fit by 

r(g) = e~ k ‘ 4/iI . (46) 

The effective absorption coefficients, that provide the best fit for each band are given 
in Table 6. 


Table 6: The effective absorption coefficient for water vapor continuum absorption in Bands 4-6 
and the three sub-bands of Band 3. Units are g ^m 2 . 


Band 3a 

Band 3b 

Band 3c 

Band 4 

Band 5 

Band 6 

£ 109.6 

54.8 

27.4 

15.8 

9.4 

7.8 


4.4 Overlapping of Absorptions 

When there is more than one absorber involved in a spectral band, overlaps must be con- 
sidered. The total Planck-weighted transmittance for two absorbers, can be written as 

_ 5 B v t i(ijij(# / 47 \ 

Tt / B u du 

where the subscripts 1 and 2 denote the two absorbers. If the transmittances due to the 
individual absorbers are expressed as the sum of the band-mean transmittance, 

[ B v r{u)du ( 48 ) 

JB v du ’ 
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( 49 ) 


and the deviation from the band-mean, r' = r - r , then (47) reduces to 


TT = TiT 2 + 


j B v T[{y)T' 2 {y)dv 
J B v dv 


If the overall shapes of the absorption curves due to both absorbers are uncorrelated with, 
each other and with B u , the second term on the right-hand side of (49) can be neglected 
and the total transmittance becomes 


tt — T iT 2 . 


(50) 


Overlapping of absorption in individual bands are shown in Table 1. As shown in Chou et 
al. (1993), the multiplication approximation (50) can be applied to Bands 4, 5, and 6 for 
the overlapping of water vapor line and continuum absorption, as well as for overlapping 
the water vapor and ozone absorption in Band 5 (9.6 /xm). Overlapping in Band 3 (15 ^m) 
is discussed in the next subsection. 


4.5 Special Treatment of the 15 fim Band 

The 15 /im band poses additional difficulties for two reasons. First, the water vapor line 
absorption and the continuum absorption are highly correlated in Band 3. As can be seen 
in the bottom two panels of Figure 3, the absorption coefficient increases with decreasing 
wavenumber by a factor of about 100 for water vapor line absorption and of about 10 for 
continuum absorption. Thus, the approximation (50) cannot be applied directly to the 
entire band. Second, the C0 2 absorption coefficients differ by several orders of magnitude 
between the band center and the wings (see the top panel of Figure 3). Rather than trying to 
parameterize the correlation effect or the variations in C0 2 absorption, we simply divide the 
band into three sub-bands (see Figure 3 and Table 1) and then combine the parameterized 
transmittances of the sub-bands into a single band transmittance. This transmittance is 
then used in the usual way to solve the transfer equations for the entire band. 


The H 2 0 Transmittance 

Water vapor line absorption in Band 3 is always done with the Jt-distribution method (Ta- 
ble 1). After dividing the band into three sub-bands, we apply the ^-distribution method 
to each sub-band as just described for the full bands. Within each of the three sub-bands 
the transmittances due to line and continuum water vapor absorption are sufficiently un- 
correlated that we can overlap them using the multiplication approximation (50). Letting 
ft an d rp be the line and continuum transmittances for sub-band i, we obtain the total 
water vapor transmittance in Band 3 by overlapping them and taking the Planck-weighted 
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Wavenumber (cm’ 1 ) 


Figure 3: Schematic of the absorption coefficient for C0 2 and for water vapor line and continuum 
absorption in the 15 /mi region (Band 3). The subdivision of the band into three sub-bands is made 
to accurately account for the large variations of the absorption coefficients. 


average: 

,wv_y c (51) 

~h B{ ?°y 

where B t {T 0 ) is the Planck flux integrated over sub-band i, and B(T 0 ) = Ei=l B i( T o ) 1S the 
total Planck flux for Band 3. 

Substituting for the line transmittance from (34) and for the continuum transmittance from 
(46) we have: 

r wv = £ e- fc ^ £ ^ fcn, 7 P (As)n,f^ 

i—l n=l ^ 

= £ e “ W £ e- k ^(A g)n,i. (52) 

1=1 n - 1 

Since we use the same scaling parameters (Table 3 and (43)) for the three sub-bands, the 
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scaled water vapor amounts, q and q are independent of sub-band. In addition, since we 
use a single set of k n , the exponentials e ~ kn */** in the inner sum of (52) only need to be 
evaluated once. 

Values for (A^) nii are given in Table 7. The effective continuum absorption coefficients, i fef, 
are given in Table 6. Note that the coefficients for Sub-bands 3a and 3b are integer multiples 
of the coefficient for Sub-band 3a; thus only one exponentiation needs to be performed to 
evaluate (46) for the three sub-bands. 


The CO 2 Transmittance 

Unlike the water vapor line transmittance, the C0 2 transmittance may be computed either 
using the ^-distribution method or by table look-up. When the table look-up is used, a 
single transmittance is computed directly for the entire band. When the ^-distribution 
method is used, we separate Sub-band 3b (the center region) from sub-bands 3a and 3c 
(the wings). The optical properties of Sub-bands 3a and 3c are very similar, and so we have 
combined them by using the same scaling parameters (a, /?, p T and m in Table 3) and the 
same set of absorption coefficients (given by k x and 77 in Table 7). The band-averaged C0 2 


Table 7: Parameters for computing transmittances in the three sub-bands of Band 3. The factors 
Ag n appear in (52) and (53). Units of k are g -1 cm 2 for water vapor and (cm-atm)J^ p for C0 2 . 


Water Vapor 


co 2 



Band 3a 

Band 3b 

Band 3c 

Wings 

Center 

b. 

n 

1.33e-2 

1.33e-2 

1.33e-2 

2.66E-5 

2.66E-3 

JL 

8 

8 

8 

8 

8 


.0000 

.0923 

.1782 

.1395 

.0766 

A<7 2 

.1083 

.1675 

.0593 

.1407 

.1372 


.1581 

.0923 

.0215 

.1549 

.1189 

A £4 

.0455 

.0187 

.0068 

.1357 

.0335 

A ^ h 

.0274 

.0178 

.0022 

.1820 

.0169 

A^6 

.0041 

.0000 

.0000 

.0220 

.0059 
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transmittance for Band 3 is thus computed from 

r C0, = £ e -( fc c)n9=/P(A ^ c ) n + jr e -(*->«fc/*(A^) nl (53) 

n=l n=l 

where subscript c denotes the band “center” (Band 3b) and subscript w the wings (Bands 
3a and 3c). Values of (A ~g) n are given in Table 7. 

Finally, the total flux transmittance in the 15 pm band is computed from 

t = t wv t co 3 ( 54 ) 


5 Vertical Discretization 


To approximate the vertical integrals in (20a, b), the atmosphere is divided into L layers 
numbered as shown in Figure 4. The upward and downward fluxes at level l for the ith 
band can be computed from 


F h 

L+l 

= i=i.. 

l'=l ’ 2 

+ 1, 

(55a) 

F h 

= E® .. 4.1 1<(''+ 1 .')-t'( 1 '.1)]. ' = 2.- 

' t,* "r a 

V=1 2 

..,L + 1, 

(55b) 


where the integer subscripts denote atmospheric levels and half-integer subscripts the layers 
between them, at which the atmospheric state variables, such as temperature and specific 
masses are defined. The quantity B i ll+ l is the spectrally integrated Planck flux of the layer 
between /' and /' + 1, and is the Planck-weighted flux transmittance between l and 

l'. The earth’s surface (level L + 1) is treated as if it were an atmospheric layer filled with 
black clouds, i.e. , B i — B{[T S ) and t*(L 4- 2, /) = 0 for I — 1. 


Because t*(1',1) is a symmetric matrix, calculations of the transmittance can be reduced 
in half by evaluating this matrix before performing the vertical integrations. If we did 
this, then in a vectorized computer code, where fluxes are computed simultaneously for m 
atmospheric columns, the l) matrix would require storage for f(L + 1)L floating point 

numbers. This could reach 10 6 storage locations. To avoid this large storage requirement, 
(55a, b) are rearranged so that the calculation can be done with only m storage locations. 


By expanding (55a, b) and rearranging terms, the upward and downward fluxes can be 
written as 


F 






L + 1 

+ E ' 


-:(i\ o 




- B W-l 


1 = 


(56a) 
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1 


2 


1-1 



Figure 4: The vertical grid and placement of various quantities for an atmosphere consisting of 
L layers. Quantities defined at the layers, such as the Planck flux B, are denoted by half-integer 
subscripts, and quantities defined at the levels separating them, such as the upward and downward 
fluxes F, by integer subscripts. The transmittance shown is for a multi-layer region bounded by 
levels l and V . Note that the surface is treated as a fictitious layer at L -f |. 
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Z'-l 

‘■‘"2 1=1 






l' = 2,...,L+1, 


(56b) 


with the boundary conditions: 

= HT.), ( 57 ») 

F} tl = 0. (57b) 

As opposed to (55a, b), only one transmittance appears in each term under the summation 
in (56a,b). The sums can thus be built-up term by term. When the transmittance of a layer 
bounded by l and V is computed, the upward flux at the top of the layer and the downward 
flux at the bottom of the layer are immediately updated; in this way there is no need to 
store the transmittance matrix, The storage for the entire routine scales then like L, rather 
than L 2 — a very significant advantage for models with high vertical resolution. 


In practice, we do not store both the upward and downward fluxes, but keep only the total 
net downward flux summed over the bands: 


i? et = - *&)■ (58) 

t=l 

This is computed for both "all-sky” and "clear-sky” conditions. The clear-sky flux is ob- 
tained from (56a, b) by replacing with the gaseous transmittance, r(l',l). 

The parameterization also computes the rate of change of the L + 1 all-sky fluxes with 
respect to the surface temperature, T t . This is done because in the general circulation 
model the IR parameterization is called relatively infrequently (currently every 3 hours), 
while the boundary layer and land surface parameterizations use time steps of a few minutes. 
To maintain consistency between the upward radiative fluxes aloft and that at the surface 
when the surface temperature changes significantly between call to the IR parameterization, 
all fluxes are linearized about the surface temperature at the beginning of the radiation 
interval, and radiative heating rates are recomputed based on this linearization every few 
minutes, at each time step. The partial derivative of the net downward flux with respect to 
the surface temperature is: 


dFp et 

dT, 


dF? 

' dT a 


= -X>*(L + M)- 


dB; 




dT a 


I = 1 1 , 


(59) 


3B * 

where r*(L + 1,1 + 1) = 1, B iL+ 1 = B t {T a ), and -fc*- is obtained by differentiating 

( 22 ). 
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6 Comparisons With Line-by-Line Calculations 


In this section, we compare fluxes and cooling rates computed from the transmittance 
parameterizations with detailed line-by-line calculations. Our line-by-line calculations of 
the absorption coefficient use the Air Force Geophysical Laboratory (AFGL) 1992 edition 
of the molecular line parameters (Rothman et al. 1987). The line profile is assumed to follow 
the Voigt function. The absorption coefficient is taken to be zero at wavenumbers > 10 
cm -1 from the line center; this is equivalent to a line-cutoff of 10 cm -1 . Since the Doppler 
line-width increases linearly with wavenumber, the spectral resolution can be increased at 
the higher wavenumbers. We use spectral intervals of 0.0005 cm” 1 for Bands 1 and 2, 
0.001 cm -1 for Band 3, and 0.002 cm” 1 for Bands 4-8. This resolution adequately resolves 
individual absorption lines. As in Ridgway et al. (1991), the flux transmittance given by 
(5) is evaluated using table look-up. 

Fluxes and cooling rates are computed for a midlatitude summer (MLS) atmosphere and a 
sub-arctic winter (SAW) atmosphere taken from McClatchey et al. (1972). Following the 
specifications of the Intercomparison of Radiation Codes in Climate Models (ICRCCM) 
(Ellingson et al. 1991), CO 2 concentration is fixed at 300 ppmv, and the specific humidity 
above the tropopause is set to 4xl0~ 6 . The atmosphere is divided into 75 layers with 5p ^ 
25 mb at pressures greater than 100 mb and 6 log 10 p = 0.15 at pressures less than 100 mb. 
In interpolating from the levels provided in McClatchey et al. (1972) to the levels used for 
the calculations, we assume the temperature, ozone concentration, and the logarithm of 
specific humidity vary linearly in the logarithm of pressure. 

The total cooling rates computed using the line-by-line method (solid lines), the HIGH 
option of the parameterization (solid circles), and the LOW option of the parameterization 
(dashed lines), are shown in Figure 5. The top panels are for the midlatitude summer 
atmosphere and the bottom panels are for the sub-arctic winter atmosphere. Similar results 
for the individual bands are shown in Figures 8 - 15 in Appendix A. 

The cooling rate error for the total IR spectrum (Figure 5) is w 0.4 C day” 1 in the mid- 
dle atmosphere and ^ 0.2 C day” 1 in the troposphere and lower stratosphere. For the 
LOW option, the cooling rate is computed accurately in the lower atmosphere, at pressures 
greater than ^ 20 mb; above the 20 mb pressure level, the cooling rate computed using the 
parameterization decreases rapidly to nearly zero. This is due to the pressure scaling of 
the absorption, which greatly underestimates the cooling in low pressure regions. It can be 
seen in the figures in the appendix that for the HIGH option, the cooling rate is computed 
accurately not only for the total cooling but also for individual bands. The maximum error 
of « 0.4 C day” 1 occurs in Band 3 (Figure 10) at pressures less than 1 mb. This error is 
small when compared with the maximum cooling of 10-12 C day” 1 due to CO 2 at these 
levels. Cooling rates in bands dominated by water vapor and ozone absorption are very 
accurately computed, with errors of less than 0.2 C day” 1 . 
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0-3000 cm 


Mid-Latitude Summer 




Figure 5: Cooling rates computed using the line-by-line method (solid lines), the “LOW” option 
(dashed lines), and the “HIGH” option (solid circles) for the total IR spectral region (0-3000 cm' 1 ). 
Cooling is due to water vapor molecular line and continuum absorption, as well as C0 2 and 0 3 
absorption. 
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The downward fluxes at the surface and the upward fluxes at the top of the atmosphere 
for the entire IR spectrum are shown in Table 8. In addition to our line-by-line results and 
the results for the LOW and HIGH options of the parameterization, we give the results 
of two other line-by-line calculations that were performed at the Goddard Laboratory for 
Atmospheres (GLA) and at the Geophysical Fluid Dynamics Laboratory (GFDL) as part 
of the ICRCCM intercomparison. These are taken from Table 2 of Ridgway et al. (1991). 
These two line-by-line calculations are very similar to each other, but differ from the line- 
by-line calculation used for the current study in important details. 


Table 8: Downward fluxes at the surface, fL , and upward fluxes at the top of the atmosphere, F ^ Q p , 
computed using a line-by-line method and the HIGH and LOW options of the parameterization. Also 
shown are the fluxes from other line-by-line calculations (GLA and GFDL) taken from, Ridgway et 
al. (1991). Units are W m“ 2 . 




F 1 

•‘top 

MID-LATITUDE SUMMER 



Parameterization (LOW) 

342.33 

292.49 

Parameterization (HIGH) 

342.45 

293.07 

line-by-line 

339.93 

293.10 

line-by-line (GLA) 

343.07 

290.45 

line-by-line (GFDL) 

343.27 

290.32 

SUB-ARCTIC WINTER 



Parameterization (LOW) 

161.17 

203.58 

Parameterization (HIGH) 

160.44 

204.22 

line-by-line 

161.51 

204.39 

line-by-line (GLA) 

164.65 

203.42 

line-by-line (GFDL) 

164.44 

202.95 


Comparing the parameterized results with our own line-by-line calculation, we find that 
the parameterizations can compute the fluxes to within 1% of the line-by-line calculations 
for the midlatitude summer atmosphere and within 0.5% for the sub-arctic winter atmo- 
sphere. Tables 9 and 10 in Appendix B give the line-by-line and parameterized fluxes for 
the individual bands. Errors are generally much smaller than 1 W m~ 2 . 

As may be seen from Table 8, there is generally better agreement between the parameteri- 
zation and our own line-by-line calculation than there is between our line-by-line calculation 
and those of GLA and GFDL; these two, on the other hand, are in very close agreement with 
each other. This disagreement between our line-by-line calculations — on which the param- 
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eterizations are based— and those used for ICECCM were somewhat disturbing. Since the 
GLA calculation was done with essentially the same code as ours, but using the conditions 
prescribed by the intercomparison, we could easily isolate the differences. Figure 6 shows 
the differences for the upward and downward fluxes for the MLS and SAW atmospheres 
between the GLA and our line-by-line calculations. Differences in Bands 1, 5, and 7 are 
small. Differences in bands 2, 3, 4, and 8 are of opposite sign at the top and bottom of the 
atmosphere, with the upward flux at the top being larger in our calculations than in the 
GLA calculations; this indicated some missing absorption in our calculations m these bands. 
This effect is particularly large in Band 2. The difference in Band 6 is more perplexing. 
These differences can be explained as follows: 


• The water vapor continuum absorption in the spectral region with strong molecular 
line absorption is not adequately understood. We therefore do not include the water 
vapor continuum absorption in the spectral region with wavenumber less than 540 
cm -1 , but the water vapor continuum absorption in the spectral region 400-540 cm 
was included in the previous GLA and GFDL calculations. This discrepancy accounts 
for most of the discrepancies in Band 2. 

• The absorption due to 0 3 in the 14 pm region, which was included in the previous GLA 
and GFDL calculations but not in the present study, contributes to the discrepancy 
in Band 3. 

• The minor C0 2 bands located at 10.4 pm (in Band 4) and 9.4 pm (in Band 5) regions 
were included in the previous GLA and GFDL calculations but not in the present 
study According to the line-by-line calculations of Kratz et al. (1993), each of the 
two minor bands causes a reduction of 0.2-0.3 W m" 2 in the top-of-the-atmosphere 
flux and an increase of 0.4-0. 5 W m -2 in the surface flux. 

• The absorption due to C0 2 in the 4.3 pm band (in Band 8) was included in the 
previous GLA and GFDL calculations but not in the present study. Most of the 
discrepancies in Band 8 as shown in Figure 6 are due to the C0 2 absorption. 

• The current study uses the 1992 edition of the AFGL line parameters, while the 
previous GLA and GFDL calculations used the 1986 edition. Integrated over 10 
cm" 1 spectral regions, the molecular line intensity of the 1992 edition is significantly 
higher (with a maximun of 33%) than the 1986 edition in the 950-1500 cm" 1 spectral 
region (Bands 5 and 6). The large difference in the downward surface fluxes for the 
midlatitude summer atmosphere is due to the use of different editions of the line 
parameters. 


These comparisons between the parameterized and line-by-line results and between the 
various line-by-line calculations give some idea of the sources and magnitudes of the errors 
that remain in the parameterization of the clear-sky IR radiation. As we can see from Table 
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DOWNWARD FLUX AT SURFACE 



* -1 


-2r -I 

- 3 i 1 1 i i i i i i * 

1 2 3 4 5 6 7 8 

BAND 

Figure 6: Differences between the GLA and our line-by-line calculations for the downward flux at 
the surface and the upward flux at the top of the atmosphere for the midlatitude summer (MLS) 
and sub-arctic winter (SAW) atmospheres. Positive values indicate our fluxes are larger than the 
GLA fluxes. 
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8, errors due to the broad-band parameterization of the line-by-line results are comparable 
to the differences between the line-by-line results. The effects of neglecting the minor bands 
o fC0 2 and O3 and of the uncertanty in modeling the water vapor continuum are significant. 
We have also neglected the effects of minor infrared absorbers, most notably nitrous oxide 
(N 2 0), methane (CH 4 ), and the chlorofluorocarbons (CFCs). The total effect of neglecting 
all the' minor absorption bands is an underestimate of « 5 W m -2 in the downward flux 
at the surface and an overestimate of * 3 W m" 2 in the upward flux at the top of the 
atmosphere. We plan to include both the minor gases and the minor bands of the mam 
gases in future versions of the parameterization. 


7 Summary 


To achieve both requirements in accuracy and speed, various approaches are applied to 
computing the transmission functions in various sections of the spectrum due to different 
absorbers. There are two options for flux calculations, HIGH and LOW. For the HIGH 
option, the /.-distribution method with linear pressure-scaling is applied to those water vapor 
spectral regions where absorption is not strong and where the contribution to the middle 
atmospheric cooling is negligible. For the C0 2 and 0 3 bands, as well as the water vapor 
spectral regions with strong absorption, cooling in the middle atmosphere is significant, 
and the transmittances are derived from pre-computed tables. For the LOW option, the 
Jt-distribution method is used in all situations, except for the absorption due to ozone which 
uses the table look-up method. 

The HIGH option of the parameterization can compute accurately the cooling rate for^the 
middle and lower atmospheres from 0.01 mb to the surface. Errors are < 0.4 C day in 
cooling rates and < 1% in fluxes. The LOW option is computationally faster than the HIGH 
option. It can achieve the same accuracy as the HIGH option except m the regions with 
pressure < 20 mb. For the HIGH option of the parameterization, it takes 9.1 seconds to 
compute the fluxes at 75 levels for 1000 soundings on a single processor of the CRAY/C90 
computer. For comparison, it takes 4.6 seconds to do the same computation with the LOW 
option. Figure 7 shows each band’s contribution to the total time. Bands 4, 5, 6, and 8 
are the same for the two options. Bands 4 and 6 differ from Band 8 only in the presence of 
continuum absorption, which evidently does not contribute significantly to the total cost. 
In Bands 1, 2, and 7, the two options differ only in the treatment of water vapor line 
absorption, ’which costs nearly four times as much when done with table look-up as when 
done by the fc-distribution method. With the HIGH option, Band 3 is the most costly due 
to the use of a table look-up to compute the C0 2 transmittances and to the separation into 
three sub-bands used to compute the water vapor transmittances. For the LOW option, 
Band 5 is the most costly due to the table look-up used to compute the 0 3 transmittances. 
A LOW option without 0 3 , which should be acceptable for tropospheric models, would be 

considerably faster. 
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Figure 7: Execution times on a single Cray C90 processor. 


In addition to the off-line testing described here, the parameterizations developed at God- 
dard have been implemented and tested in long-term simulations in the climate and data 
assimilation general circulation models at Goddard. 
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Appendix A: Band-by-Band Reults 

0-340 cm' 1 


Mid-Latitude Summer 




Figure 8: Cooling rates computed using the line-by-line method (solid lines), the “LOW” option 
(dashed lines), and the “HIGH” option (solid circles) for the 0-340 cm -1 band. Cooling is due to 
water vapor molecular line absorption. Options for the transmittance parameterization are outlined 

in Table 1 
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800-980 cm' 1 


Mid-Latitude Summer 



Figure 11. Same as Figure 8, except for the 800—980 cm * band. Cooling is due to water vapor 
molecular line and continuum absorption. Note that the HIGH and LOW options are identical in 
this band. 
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11 00-1 380 cm' 1 


Mid-Latitude Summer 



Sub-Arctic Winter 



Figure 13: Same as Figure 8, except for the 1100-1380 cm 1 band. Cooling is due to water vapor 
molecular line and continuum absorption. 
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Pressure (mb) Pressure (mb 


1 





Table 9: Downward fluxes at the surface, F^ c , and upward fluxes at the top of the atmosphere, 
Fi op , for the mid-latitude summer atmosphere computed using a line-by-line method and the HIGH 
and LOW options of the parameterization. Units are W m -2 . 



Spectral Band (cm l ) 


*\ T op 

0 

- 340 




line-by-line 

50.96 

34.25 


Parameterization (HIGH) 

50.97 

34.04 


Parameterization (LOW) 

50.97 

33.92 

340 

- 540 




line-by-line 

80.72 

60.51 


Parameterization (HIGH) 

81.23 

60.01 


Parameterization (LOW) 

81.28 

60.03 

540 

- 800 




line-by-line 

105.98 

68.03 


Parameterization (HIGH) 

107.58 

68.40 


Parameterization (LOW) 

107.43 

67.74 

800 

- 980 




line-by-line 

27.97 

58.49 


Parameterization (HIGH) 

28.35 

58.50 


Parameterization (LOW) 

28.34 

58.50 

980 

- 1100 




line-by-line 

12.80 

22.28 


Parameterization (HIGH) 

12.86 

21.81 


Parameterization (LOW) 

12.86 

21.82 

1100 

- 1380 




line-by-line 

28.14 

37.18 


Parameterization (HIGH) 

27.95 

38.21 


Parameterization (LOW) 

27.95 

38.21 

1380 

- 1900 




line-by-line 

30.30 

7.27 


Parameterization (HIGH) 

30.35 

7.22 


Parameterization (LOW) 

30.33 

7.40 

1900 

- 3000 




line-by-line 

3.06 

5.09 


Parameterization (HIGH) 

3.16 

4.88 


Parameterization (LOW) 

3.16 

4.88 
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Table 10: Same as Table 9, except for the sub-arctic winter atmosphere. 


Spectral Band (cm l ) 

F * 

r *fc 

*\ T op 

0 - 340 

line-by-line 

40.39 

32.10 

Parameterization (HIGH) 

40,40 

31.93 

Parameterization (LOW) 

40.40 

31.82 

340 - 540 

line-by-line 

47.35 

52.01 

Parameterization (HIGH) 

47.30 

51.74 

Parameterization (LOW) 

48.09 

51.78 

540 - 800 

line-by-line 

53.17 

51.38 

Parameterization (HIGH) 

51.84 

51.65 

Parameterization (LOW) 

51.90 

51.06 

800 - 980 

line-by-line 

1.45 

32.85 

Parameterization (HIGH) 

1.63 

32.86 

Parameterization (LOW) 

1.63 

32.86 

980 - 1100 

line-by-line 

3.14 

10.99 

Parameterization (HIGH) 

3.23 

10.87 

Parameterization (LOW) 

3.22 

10.87 

1100 - 1380 

line-by-line 

5.49 

18.80 

Parameterization (HIGH) 

5.49 

18.98 

Parameterization (LOW) 

5.50 

18.98 

1380 - 1900 

line-by-line 

10.11 

4.90 

Parameterization (HIGH) 

10.13 

4.87 

Parameterization (LOW) 

10.02 

4.90 

1900 - 3000 

line-by-line 

0.41 

1.36 

Parameterization (HIGH) 

0.42 

1.32 

Parameterization (LOW) 

0.43 

1.32 
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Appendix B: The Code 


The following is a listing of the FORTRAN implementation of the parameterization. This 
code follows the “plug-compatible” rules of Kalnay et al. (1988), which should make it easy 
to use in climate models. The parameterization is accessed by calling subroutine IRRAD. 
All information needed by this subroutine is passed through the argument list and all 
results are returned through the argument list. Its arguments are described in detail in the 
comments at the top of the code. 


The inputs to the routine are a two-dimensional (M by N) array of soundings, with each 
sounding consisting of NP layers. The temperature, the specific humidity, the ozone mixing 
ratio, and the cloud fraction and optical thickness are specified at the NP layers, while 
the pressure is specified at the NP-1 interfaces and at the lower and upper boundaries. 
The surface temperature is also required. The CO 2 concentration is assumed constant 
throughout the atmosphere. 

Although the routine processes M by N soundings, all input and output arrays are dimen- 
sioned M by NDIM, where NDIM > N. This is done so that a large array can be stripmined 
over the outer dimension without the need to make copies. Such stripmining is useful in 
parallel applications and when trying to minimize storage, since the scratch space required 
by the routine is proportional to M*N. 

The outputs of the routine are the net downward all-sky and clear-sky fluxes and the rate of 
change of the aU-sky flux with respect to surface temperature. The routine also returns the 
upward flux at the surface, i.e., the approximate form of crTf used in the parameterization. 
This allows the calling program to maintain energetic consistency by using this value in the 
surface energy budget. 
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7.1 Subroutine IRRAD 


********************* CLIRAD IR Date: October, 1994 **************** 
* 

subroutine irrad (m,n,ndim,np,taucl,ccld,pl,ta,wa,oa,co2,ts, 

* high,f lx,flc,df dts,st4) 

* 


*****************************4 C ******* ++ * + *,i t j* ++J | t++++++J|t+1(lJt[J | {1Jtl j {1J{J | {) i, ##+J j tJ((1 j t+ 

* 

* This routine computes ir fluxes due to water vapor, co2, and o3. 

* Clouds in different layers are assumed randomly overlapped. 

* 

* This is a vectorized code. Xt computes fluxes simultaneously for 

* (m x n) soundings, which is a subset of (m x ndim) soundings. 

* In a global climate model, m and ndim correspond to the numbers of 

* grid boxes in the zonal and meridional directions, respectively. 

* 

* Detailed description of the radiation routine is given in 

* Chou and Suarez (1994). 

* 

* There are two options for computing cooling rate profiles. 

* 

* if high « .true., transmission functions in the co2, o3, and the 

* three water vapor bands with strong absorption are computed using 

* table look-up. cooling rates are computed accurately from the 

* surface up to 0.01 mb. 

* if high » .false., transmission functions are computed using the 

* k-distribution method with linear pressure scaling, cooling rates 

* are not calculated accurately for pressures less than 20 mb. 

* the computation is faster with high=. false, than with high^.true. 

* 

* The IR spectrum is divided into eight bands: 


wavenumber (/cm) absorber 


method 


* 

1 

0 

- 340 

h2o 

K/T 

* 

2 

340 

- 540 

h2o 

K/T 

* 

3 

540 

- 800 

h2o ,cont ,co2 

K.S.K/T 

* 

4 

800 

- 980 

h2o ,cont 

K.S 

* 

5 

980 

- 1100 

h2o ,cont ,o3 

K.S.T 

* 

6 

1100 

- 1380 

h2o,cont 

K.S 

* 

7 

1380 

- 1900 

h2o 

K/T 

* 

8 

1900 

- 3000 

h2o 

K 


* 


* Note : "h2o" for h2o line absorption 

* "cont" for h2o continuum absorption 

* "K" for k-distribution method 

* "S" for one-parameter temperature scaling 

* "T" for table look-up 

* 


♦ The 15 micrometer region (540-800/cm) is further divided into 

* 3 sub -bands : 
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* 

* 

* 

♦ 

* 

* 


subbnad wavenumber (/cm) 


540 - 620 
620 - 720 
720 - 800 


*— 

— Input parameters 

units 

size 

♦ 

♦ 

number of soundings in zonal direction (m) 

n/d 

1 

* 

number of soundings in meridional direction (n) 

n/d 

1 

* 

* 

maximum number of soundings in 

meridional direction (ndim) 

n/d 

1 

* 

number of atmospheric layers (np) 

n/d 

1 

* 

cloud optical thickness (taucl) 

n/d 

m*ndim*np 

* 

cloud cover (ccld) 

fraction 

m*ndim*np 

* 

level pressure (pi) 

mb 

m*ndim*(np+l) 

* 

layer temperature (ta) 

k 

m*ndim*np 

* 

layer specific humidity (wa) 

g/g 

m*ndim*np 

* 

layer ozone mixing ratio by mass (oa) 

g/g 

m*ndim*np 

* 

surface temperature (ts) 

k 

m*ndim 

* 

co2 mixing ratio by volumn (co2) 

ppmv 

1 

* 

high 


1 


pre-computed tables used in table look-up for transmittance calculations. 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 


cl , c2 , c3 : for co2 (band 3) 
ol , o2, o3 : for o3 (band 5) 
hi 1 ,hl2 ,hl3 : for h2o (band 1) 
h21,h22,h23: for h2o (band 2) 
h71,h72,h73: for h2o (band 7) 

output parameters 

net downward flux, all-sky (fix) 
net downward flux, clear-sky (flc) 
sensitivity of net downward flux 
to surface temperature (dfdts) 
emission by the surface (st4) 


w/m**2 

w/m**2 


m*ndim* (np+1) 
m*ndim* (np+1) 


w/m**2/k m*ndim*(np+i) 
w/m**2 m*ndim 


Notes : 

( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


(7) 

( 8 ) 

(9) 


Water vapor continuum absorption is included in 540-1380 /cm. 
Scattering by clouds is not included. 

Clouds are assumed "gray 1 ' bodies. 

The diffuse cloud transmission is computed to be exp(-l . 66*taucl) . 

If there are no clouds, f lx=f lc . 

plevel(l) is the pressure at the top of the model atmosphere, and 
plevel(np+l) is the surface pressure. 

Downward flux is positive, and upward flux is negative, 
dfdts is always negative because upward flux is defined as negative. 
For questions and coding errors, please contact with Ming-Dah Chou, 
Code 913, NASA/Goddard Space Flight Center, Greenbelt, MD 20771. 


44 


Phone: 301-286-4012, Fax: 301-286-1759 
e-mail : chouOclimate.gsfc.nasa.gov 


* 

* 

* 

c parameters defining the size of the pre-computed tables for transmittance 
c calculations using table look-up. 
c 

c "m" is the number of intervals in pressure 
c "no" is the number of intervals in o3 amount 

c "nc" is the number of intervals in co2 amount 

c "nh" is the number of intervals in h2o amount 

c "nt" is the number of copies to be made from the pre-computed 

c transmittance tables to reduce "memory-bank conflict" 

c in parallel machines and, hence, enhancing the speed of 

c computations using table look-up. 

c If such advantage does not exist, "nt" can be set to 1. 

c***** *************************************** **************************** 41 ** 
implicit none 

integer nx,no,nc,nh,nt ,m,n,ndim,np 
integer i, j ,k,ip,iw,it,ib,ik,iq,isb,kl,k2 
parameter (nx=26 ,no=2i ,nc=24 ,nh=3i ,nt=7) 

c input parameters 

real*8 co2 

real*8 taucl (m,ndim,np) ,ccld(m,ndim,np) ,pl(m,ndim,np+l) , 

* ta(m,ndim,np) , wa(m,ndim,np) ,oa(m,ndim ,np) ,ts(m,ndim) 

logical high 

c output parameters 

real*8 f lx (m,ndim,np+l) ,f lc (m,ndim,np+l) ,df dts (m,ndim,np+l) , 

* st4(m,ndim) 

c static data 

real*8 cb(5,8) 

c temporary arrays 

real*8 pa(m,n,np) ,dt (m,n,np) 

real*8 sh2o (m,n ,np+l) ,swpre(m,n,np+l) , swtem(m,n,np+l) 

real*8 sco3 (m,n,np+l) ,scopre(m,n,np+l) , scotem(m,n,np+l) 

real*8 dh2o(m,n,np) ,dcont (m,n,np) ,dco2(m,n,np) ,do3(m,n,np) 

real*8 th2o(m,n,6) ,tcon(m,n,3) , tco2(m,n f 6 , 2) 

real*8 h2oexp (m,n,np , 6) , conexp(m,n,np, 3) , co2exp(m,n,np, 6, 2) 

real*8 clr (m,n,0 :np+l) ,fclr (m,n) 

real*8 blayer (m,n, 0 : np+1) ,dbs(m,n) 

real*8 trant(m,n) 

logical oznbnd 
logical co2bnd 
logical h2otbl 
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logical conbnd 


real*8 cl (nx.nc ,nt) ,c2 (nx,nc,nt) ,c3 (nx,nc,nt) 
real*8 ol (nx ,no ,nt) ,o2 (nx ,no ,nt) ,o3 (nx.no.nt) 
real*8 hll (nx,nh,nt) ,hl2(nx,nh,nt) ,h!3 (nx ,nh,nt) 
real*8 h21 (nx,nh,nt) ,h22(nx,nh,nt) ,h23 (nx ,nh,nt) 
real*8 h71 (nx,nh,nt) ,h72(nx,nh,nt) ,h73(nx,nh,nt) 

real *8 dp,xx,trantc 

— the following coefficients (table 2 of chou and suarez, 1995) are for 
computing spectrally integtrated planck fluxes of the 8 bands using 

( 22 ) 

data cb/ 

1 -2.6844e-l,-8.8994e-2. 1 . S676e-3 , -2 . 9349e-6 , 2.2233e-9, 

2 3 . 7315e+i ,-7.4758e-l, 4. 6151e-3 ,-6 .3260e-6, 3.5647e-9, 

3 3.7187e+i,-3.9085e-l,-6. 1072e-4, 1 .4534e-5 ,-l .6863e-8, 

4 -4 . 1928e+l , 1 . 0027e+0,-8 . 5789e-3, 2.91996-6,-2.56546-8, 

5 -4 . 9163e+l , 9 . 8457e-l ,-7 . 0968e-3, 2.0478e-5,-l .5514e-8, 

6 -1.0345e+2, 1 . 8636e+0, -1 . 1753e-2 , 2 . 7864e-5,-i . 1998e-8, 

7 -6 . 9233e+0,-i . 5878e-l , 3 . 9160e-3,-2 .4496e-5, 4.9301e-8, 

8 1 . 1483e+2,-2.2376e+0, 1 . 6394e-2, -5 .3672e-5, 6.6456e-8/ 

— copy tables to enhance the speed of co2 (band 3) , o3 (bands) , 
and h2o (bands 1, 2, and 7 only) transmission calculations 
using table look-up. 


logical first 
data first /.true./ 

include "h2o.tran3" 
include "co2 . tran3" 
include "o3.tran3" 

save cl,c2,c3,ol,o2,o3 

save hll ,hl2 ,hl3 ,h21 ,h22 ,h23 ,h71 ,h72 ,h73 
if (first) then 

-tables co2 and h2o are only used with 'high' option 

if (high) then 

do iw=l,nh 
do ip=l,nx 

hll (ip, iw , i)=l .0-hll(ip,iw , 1) 
h21 (ip , iw , 1) =1 .0-h21 (ip, iw , 1) 
h71(ip,iw,l)=1.0-h71(ip,iw,l) 
enddo 
enddo 


do iw=i,nc 
do ip=i,nx 

cl (ip, iw, l)=i.O-cl(ip,iw,l) 
enddo 
enddo 

-tables are replicated to avoid memory bank conflicts 

do it=2,nt 
do iw=l,nc 
do ip=i,nx 

cl (ip, iw, it)= cl(ip,iw,l) 
c2 (ip,iw,it)= c2(ip,iw,l) 
c3 (ip, iw, it)= c3 (ip, iw, 1) 
enddo 
enddo 

do iw=l ,nh 
do ip=l,nx 

hll (ip, iw, it)=hll(ip,iw,l) 
hi 2 (ip, iw, it)=hl2(ip, iw, 1) 
hl3(ip,iw,it)=hl3(ip,iw, 1) 
h21 (ip, iw, it ) =h2 1 ( ip , iw , 1 ) 
h22 (ip, iw, it)=h22(ip,iw, 1) 
h23(ip, iw, it)=h23(ip, iw, 1) 
h71 (ip, iw, it)=h71(ip, iw, 1) 
h72(ip, iw, it)=h72(ip, iw, 1) 
h73(ip, iw, it)=h73(ip, iw, 1) 
enddo 
enddo 
enddo 

endif 

always use table look-up for ozone transmittance 

do iw*l,no 
do ip=l,nx 

ol (ip, iw, 1)=1 . 0-ol(ip, iw, 1) 
enddo 
enddo 

do it=2,nt 
do iw=l,no 
do ip=l ,ni 

01 (ip,iw,it)= ol(ip,iw,l) 

02 (ip,iw,it)= o2(ip,iw,l) 

03 (ip, iw,it)~ o3(ip, iw, 1) 
enddo 

enddo 

enddo 

first=. false . 
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endif 


compute layer pressure (pa) and layer temperature minus 250K (dt) 

do k= 1 , np 
do j=l # n 
do i=l,m 

pa(i, j ,k)=0.5*(pl(i, j ,k)+pl(i, j ,k+l)) 
dt(i, j ,k)=ta(i, j ,k) -250.0 
enddo 
enddo 
enddo 

compute layer absorber amount 
dh2o : water vapor amount (g/cm**2) 

dcont : scaled water vapor amount for continuum absorption (g/cm**2) 

dco2 : co2 amount (cm-atm)stp 

do3 : o3 amount (cm-atm)stp 

the factor 1.02 is equal to 1000/980 

factors 789 and 476 are for unit conversion 

the factor 0.001618 is equal to 1 . 02/ (. 622*1013 .26) 

the factor 6.081 is equal to 1800/296 

do k=I ,np 
do j=i,n 
do i=l,m 

dp = pi (i , j » k+l)~pl(i i j ,k) 

dh2o (i,j,k) = i.02*wa(i, j ,k)*dp+l . e-10 

dco2 (i,j,k) " 789. *co2*dp+l. e-10 

do3 (i, j ,k) = 476 . 0*oa(i , j ,k) *dp+l .e-10 

-compute scaled water vapor amount for h2o continuum absorption 
following eq. (43). 

rx=pa(i , j , k) *0 . 001618*wa(i, j ,k) *wa(i , j ,k) *dp 
dcont (i , j ,k) = xx*exp(1800./ta(i, j ,k) -6. 081)+1 . e-10 

-compute effective cloud-free fraction, clr, for each layer, 
the cloud diffuse transmittance is approximated by using a 
diffusivity factor of 1.66. 

clr (i, j ,k)=i . 0- (ccld(i , j ,k) * (1 . -erp(-l . 66*taucl(i, j ,k)))) 

enddo 

enddo 

enddo 

—compute column-integrated h2o amoumt, h2o-weighted pressure 
and temperature, it follows eqs . (37) and (38). 


if (high) then 


call column (m,n,np, pa, dt ,dh2o ,sh2o , swpre , swtem) 
endif 

c the surface (with an index np+1) is treated as a layer filled with 

c black clouds. 

do j-i ,n 
do i=l,m 

clr (i, j ,0) = 1.0 

clr (i, j ,np+l) =0.0 
st4(i,j) = 0.0 

enddo 
enddo 

c initialize fluxes 

do k=l,np+l 
do j=l,n 
do i=l,m 

f lx(i, j ,k) =0.0 

f lc(i, j ,k) =0.0 

dfdts (i, j ,k)= 0 . 0 
enddo 
enddo 
enddo 

c integration over spectral bands 

do 1000 ib=l , 8 

c if h2otbl, compute h2o (line) transmittance using table look-up. 

c if conbnd, compute h2o (continuum) transmittance in bands 3, 4, 5 and 6. 

c if co2bnd, compute co2 transmittance in band 3. 

c if oznbnd, compute o3 transmittance in band 5. 

h2otbl=high . and . ( ib . eq . 1 . or . ib . eq . 2 . or . ib . eq . 7) 
conbnd=ib . ge . 3 . and . ib . le . 6 
co2bnd=ib . eq . 3 
oznbnd=ib . eq. 5 

c blayer is the spectrally integrated planck flux of the mean layer 

c temperature derived from eq. (22) 

do k=l,np 
do j=l,n 
do i=l,m 

blayer(i,j,k)=ta(i,j,k)*(ta(i,j,k)*(ta(i,j ,k) 

* *(ta(i, j ,k)*cb(5,ib)+cb(4,ib))+cb(3,ib)) 

* +cb(2,ib))+cb(l, ib) 
enddo 

enddo 
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enddo 


the earth's surface, with an index "np+l”, is treated as a layer 

do j=i,n 
do i=l,m 

blayer(i, j ,0) =0.0 

blayer (i, j ,np+l)=ts(i, j)* (ts(i, j)*(ts(i, j) 

*(ts(i, j) *cb ( 5 , ib) +cb (4 , ib) ) +cb (3 , ib) ) 
+cb(2,ib))+cb(l ,ib) 

dbs is the derivative of the surface planck flux with respect to 
surface temperature (eq. 59) . 

dbs(i, j)=ts(i, j)*(ts(i, j)*(ts(i, j) 

*4 . *cb ( 5 , ib) +3 . *cb (4 , ib) ) +2 . *cb (3 , ib) ) +cb (2 , ib) 


enddo 

enddo 

compute column-* integrated absorber amoumt , absorber-weighted 
pressure and temperature for co2 (band 3) and o3 (band 5) . 
it follows eqs . (37) and (38). 

this is in the band loop to save storage 

if( high .and. co2bnd) then 

call column (m,n,np, pa, dt ,dco2 ,sco3 ,scopre,scotem) 

endif 

if(oznbnd) then 

call column (m,n,np, pa, dt ,do3 ,sco3,scopre,scotem) 
endif 

compute the exponential terms (eq. 32) at each layer for 
water vapor line absorption when k-distribution is used 

if ( .not. h2otbl) then 

call h2oexps(ib,m,n,np,dh2o,pa,dt,h2oexp) 

endif 


compute the exponential terms (eq. 46) at each layer for 
water vapor continuum absorption 

if ( conbnd) then 

call conexps (ib,m,n ,np,dcont , conexp) 
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■compute the exponential terms (eq. 32) at each layer for 
co2 absorption 

if( .not. high .and. co2bnd) then 

call co2exps(m,n,np,dco2,pa,dt ,co2exp) 
endif 

compute transmittances for regions between levels kl and k2 
and update the fluxes at the two levels. 

do 2000 kl=i,np 

initialize fclr, th2o, tcon, and tco2 

do j=l,n 
do i=l,m 
f clr(i , j)-l . 0 
enddo 
enddo 

for h2o line absorption 

if (.not. h2otbl) then 
do ik=l ,6 
do j=l,n 
do i-l,m 
th2o(i , j ,ik)=l . 0 
enddo 
enddo 
enddo 
endif 

for h2o continuum absorption 

if (conbnd) then 
do iq= 1 , 3 
do j = l ,n 
do i-i,m 
tcon(i, j , iq) =1 . 0 
enddo 
enddo 
enddo 
endif 

for co2 absorption when using k-distribution method, 
band 3 is divided into 3 sub-bands, but sub-bands 3a and 3c 
are combined in computing the co2 transmittance. 


if (.not. high .and. co2bnd) then 
do isb=l,2 
do ik=i,6 
do j^i.n 
do i=i,m 

tco2 (i, j , ik, isb) =1 . 0 
enddo 
enddo 
enddo 
enddo 
endif 

c loop over the bottom level of the region (k2) 

do 3000 k2=kl+l ,np+l 

do j«i,n 
do i=i,m 
trant (i, j)=i . 0 
enddo 
enddo 

if (h2otbl) then 

c compute water vapor transmittance using table look-up 


if (ib.eq.l ) then 

call tablup(kl,k2 l m,n,np l nx,nh,nt , sh2o , swpre ,swtem, 

-8. 0,-2. 0,0. 3, 0.2, hi 1 ,hl2 ,h!3 , trant) 

endif 

if (ib.eq.2 ) then 

call tablup(ki,k2,m,n.np,nx,nh,nt ,sh2o, swpre. swtem, 

-8. 0,-2. 0,0. 3, 0.2 ,h21 ,h22,h23 .trant) 

endif 

if (ib.eq.7 ) then 

call tablup(kl ,k2 ,m,n,np,nr,nh,nt , sh2o .swpre, swtem, 

-8. 0,-2. 0,0. 3, 0.2 ,h71 ,h72,h73, trant) 

endif 

else 

compute water vapor transmittance using k-distribution. 

call wvkdis ( ib ,m ,n ,np ,k2-l ,h2oexp , conerp , th2o , icon, trant) 


endif 
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if(co2bnd) then 


if( high ) then 

c compute co2 transmittance using table look-up method 

call tablup(kl,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem, 

* -4. 0,-2. 0,0. 3, 0.2, cl ,c2 ,c3,trant) 

else 

c compute co2 transmittance using k-distribution method 

call co2kdis(m,n,np,k2-i ,co2exp, tco2 , trant) 
endif 
endif 

c compute o3 transmittance using table look-up 

if (oznbnd) then 

call tablup(kl ,k2 ,m,n,np,nx,no ,nt ,sco3, scopre.scotem, 

* -6. 0,-2. 0,0. 3,0. 2, ol ,o2 ,o3, trant) 

endif 

c fclr is the clear line-of -sight between levels kl and k2. 

c in computing fclr, clouds are assumed randomly overlapped 
c using eq. (10) . 

do j=l,n 
do i=i,m 

f clr (i, j) = f clr(i, j)*clr(i, j ,k2-l) 
enddo 
enddo 

c compute upward and downward fluxes 

c add "boundary" terms to the net downward flux. 

c these are the first terms on the right -hand- side of 
c eqs . (56a) and (56b). 

c downward fluxes are positive. 

if (k2 .eq. kl+1) then 
do j=l,n 
do i=i,m 

f lc(i, j ,kl)=f lc(i, j ,kl)-blayer (i, j ,kl) 
f lc(i, j ,k2)=f lc (i, j ,k2)+blayer (i, j ,kl) 
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flx(i, j ,kl)=flx(i, j.kl) -blayer (i, j.kl) 
flx(i, j ,k2)=flx(i, j ,k2)+blayer (i, j ,kl) 

enddo 

enddo 

endif 

B add ilux components involving the four layers above and below 

c the levels kl and k2. it follows eqs. (66a) and (56b). 

do j=i»n 
do i=i,m 

xx=trant (i, j) * (blayer (i , j ,k2-l) -blayer (i , j »k2)) 

f lc (i» j ,kl) =flc(i, j ,kl)+xx 

flx(i, j ,kl) =flx(i, j ,kl)+xx*fclr(i, j) 

xx=trant (i , j) * (blayer (i > j . kl-1) -blayer (i, j .kl) ) 

flc(i,j.k2) =f lc(i, j ,k2)+xx 

flx(i, j.k2) =flx(i, j ,k2)+xx*f clr (i, j) 

enddo 

enddo 

3000 continue 

c- compute the partial derivative of fluxes with respect to 

c surface temperature (eq. B9) . 

do j*l,n 
do i=i f m 

dfdts(i.j.kl) =dfdts(i, j ,kl)-dbs(i, j)*trant(i, j)*fclr(i, j) 
enddo 
enddo 

2000 continue 

c add contribution from the surface to the flux terms at the surface. 

do j=l,n 
do i^l.m 

dfdts(i, j ,np+l) =df dts (i , j ,np+l) -dbs (i , j) 
enddo 
enddo 

do j=i,n 
do i=l ,m 

f lx (i, j,np+l)=flx(i, j.np+l) -blayer (i,j,np+l) 
f lc (i, j ,np+l)=flc(i, j ,np+l)-blayer(i, j ,np+l) 
st4(i.j)=st4(i,j)-blayer(i, j.np+l) 
enddo 
enddo 
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1000 continue 


return 

end 



7.2 Subroutine COLUMN 


subroutine column (m, n,np, pa, dt .sabsO, sabs, spre .stem) 
c************************************************************************** 

c compute column-integrated (from top of the model atmosphere) 

c absorber amount (sabs) , absorber-weighted pressure (spre) and 

c temperature (stem) . 

c computations of spre and stem follows 

c eqs . (37) and (38) of Chou and Suarez (1994). 

c 

c input parameters 

c number of soundings in zonal direction (m) 
c number of soundings in meridional direction (n) 
c number of atmospheric layers (np) 
c layer pressure (pa) 
c layer temperature minus 250K (dt) 
c layer absorber amount (sabsO) 
c 

c output parameters 

c column- integrated absorber amount (sabs) 
c column absorber-weighted pressure (spre) 
c column absorber- weighted temperature (stem) 
c 

c units of pa and dt are mb and k, respectively. 

c units of sabs are g/cm**2 for water vapor and (cn ” atm)stp ^^^ + ^^° 3 

implicit none 
integer m,n i np,i,j,k 

c input parameters 

real*8 pa(m,n,np) ,dt(m,n,np) , sabsO (m.n.np) 

c output parameters 

real*8 sabs(m,n,np+l) ,spre(m,n,np+l) ,stem(m,n,np+l) 

do j=l,n 
do i=l,m 
sabs (i, j » 1)=0.0 

spre(i, j,l)=0.0 
stem(i, j , 1)*0. 0 
enddo 
enddo 

do k= 1 » np 
do j*i,n 
do i=i,m 

sabs(i,j.k+l)=sabs(i, j ,k)+sabsO(i. j.k) 
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spre ( i , j , k+ 1 ) =spre (i , j , k) +pa ( i , j , k) *sabsO ( i , j , k) 
stem(i, j ,k+l)=stem(i, j,k)+dt(i, j ,k)*sabsO(i, j ,k) 
enddo 
enddo 
enddo 


return 

end 



7.3 Subroutine TABLUP 


c ********************************************************************** 
subrout ine t ablup ( k 1 , k2 , m , n , np , nx , nh , nt , s abs , spre , st em , w 1 , p 1 , 

* dwe,dpe,coef 1 ,coef 2,coef 3, tran) 

c ********************* ******************** ***************************** 
c coinput© water vapor, co2, and o3 transmit tances between levels kl and k2 
c using table look-up for m x n soundings, 
c 

c Calculations follow Eq. (40) of Chou and Suarez (1994) 
c 

c input — — — 

c indices for pressure levels (kl and k2) 
c number of grid intervals in zonal direction (m) 

c number of grid intervals in meridional direction (n) 

c number of atmospheric layers (np) 

c number of pressure intervals in the table (nx) 

c number of absorber amount intervals in the table (nh) 

c number of tables copied (nt) 

c column- integrated absorber amount (sabs) 
c column absorber amount -weighted pressure (spre) 
c column absorber amount -weighted temperature (stem) 
c first value of absorber amount (loglO) in the table (wl) 
c first value of pressure (loglO) in the table (pi) 
c size of the interval of absorber amount (loglO) in the table (dwe) 

c size of the interval of pressure (loglO) in the table (dpe) 

c pre-computed coefficients (coefl, coef2, and coef3) 
c 

c updated 

c transmittance (tran) 
c 

c Note: 

c (1) units of sabs are g/cm**2 for water vapor and (cm-atm)stp for co2 
and o3 . 

c (2) units of spre and stem are, respectively, mb and K. 

c (3) there are nt identical copies of the tables (coefl, coef2, and 

c coef3) . the prupose of using the multiple copies of tables is 

c to increase the speed in parallel (vectorized) computations. 

C if such advantage does not exist, nt can be set to 1. 

c 

c ** ************************************************* ******************* 
implicit none 

integer kl ,k2 ,m,n,np,nx,nh,nt , i , j ,k 

c input parameters 

real*8 wl ,pl ,dwe ,dpe 

real*8 sabs(m,n,np+l) ,spre(m,n,np+l) ,stem(m,n,np+l) 
real*8 coef 1 (nx,nh,nt) ,coef 2(nx,nh,nt) , coef 3 (nx,nh,nt) 

c update parameter 
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real*8 tran(m,n) 


c temporary variables — 

real*8 il»i2,i3.we,pe.fB,ip,fm,pa,pb,pc,ax,ba,bb,ti,ca > cb,t2 
integer iw,ip,nn 


c +****+* ******************* ********** ********************************** 


do j=l,n 
do i=i,m 

nn=mod(i,nt)+l 

xl=sabs(i, j ,k2)-sabs(i, j ,kl) 
x2=(spre(i, j,k2)-spre(i, j,kl))/xl 
x3=(stem(i, j ,k2)-stem(i, j,ki))/xl 

we=(loglO(xl)-wl)/dwe 

pe=(logl0(x2)-pl)/dpe 

we=max (we, wl-2 . *dwe) 
pe=max(pe,pl) 

iw=int (we+1 .5) 
ip=int Cpe+1.5) 

iw=min(iw,nh-l) 
iw=max(iw, 2) 

ip=min(ip,nx-l) 
ip=max(ip, 1) 

f w=we-f loat (iw-1) 
fp=pe-float (ip-1) 

c linear interpolation in pressure 

pa = coef 1 (ip, iw-i ,nn)*(i. -fp)+coefl (ip+1 , iw-l,nn) *fp 
pb = coef 1 (ip, iv, nn)*(l .-fp)+coef 1 (ip+l t iw, nn)*fp 
pc = coef l(ip,iw+i,nn)*(l.-fp)+coefl(ip+l,iw+l,nn)*fp 

c quadratic interpolation in absorber amount 

ax = (-pa*(l.-fw)+pc*(l.+fw)) *fw*0.5 + pb* (1 . -f w*fw) 

c linear interpolation in temperature 

ba = coef 2(ip , iw , nn) *(1 . -fp)+coef 2 (ip+1 , iw, nn)*fp 
bb = coef 2 (ip , iw+1 ,nn) * (1 . -f p)+coef 2 (ip+1 , iw+i ,nn) *fp 
tl = ba*(l . -fw) + bb*f w 

ca = coef 3 (ip , iw , nn) * (1 . -fp)+coef 3 (ip+1 , iw, nn)*fp 
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cb = coet 3 (ip , iw+1 ,nn) * ( i . -lp)+coef 3 (ip+1 , iw+1 ,nn) *f p 
t2 * ca* ( 1 . -f w) + cb*fw 

update the total transmittance between levels ki and k2 
tran(i, j)= (ax + (tl+t2*x3) * x3)*tran(i, j) 


enddo 

enddo 


return 

end 


7.4 Subroutine WVKDIS 


c **************************4^******4^****%***************************** 
subrout ine wvkdis ( ib , m ,n , np , k , h2oexp , conexp , th2o , tcon, tran) 
c ********** *****************************+*********+******************** 
c compute water vapor transmittance between levels kl and k2 for 
c min soundings using the k-distribution method, 
c 

c computations follow eqs . (34), (46), (50) 
c and (52) of Chou and Suarez (1994). 
c 

c input parameters 

c spectral band (ib) 

c number of grid intervals in zonal direction (m) 
c number of grid intervals in meridional direction (n) 
c number of levels (np) 
c current level (k) 

c exponentials for line absorption (h2oexp) 
c exponentials for continuum absorption (conexp) 
c 

c updated parameters 

c transmittance between levels kl and k2 due to 
c water vapor line absorption (th2o) 
c transmittance between levels kl and k2 due to 
c water vapor continuum absorption (tcon) 
c total transmittance (tran) 
c 

C*** ********** ******* ****** ************ ******** 4c * ******************** * * 

implicit none 
integer ib,m,n,np,k, i, j 

c input parameters 

real*8 conexp (m,n,np, 3) ,h2oexp(m,n,np, 6) 

c updated parameters 

real*8 th2o(m,n,6) ,tcon(m,n,3) ,tran(m,n) 

c static data 

integer ne(8) 

real*8 fkw(6,8) ,gkw(6, 3) 

c temporary arrays 

real*8 trnth2o 


c fkw is the pi anck- weighted k-distribution function due to h2o 

c line absorption given in table 4 of Chou and Suarez (1995). 
c the k-distribution function for the third band, fkw(*,3), is not used 
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data fkw / 

2 

3 

4 

5 

6 

7 

8 


0.2747,0.2717,0.2752,0.1177,0.0352,0.0255, 

0.1521,0.3974,0.1778,0.1826,0.0374,0.0527, 

6 * 1 . 00 , 

0.4654.0.2991,0.1343,0. 0646 , 0 . 0226 ,0.0140, 
0.5543,0.2723,0.1131,0.0443,0.0160,0.0000, 
0.1846,0.2732,0.2353,0.1613,0.1146,0.0310, 
0.0740,0.1636,0.4174,0.1783,0.1101,0.0566, 
0.1437,0.2197,0.3185,0.2351,0.0647.0.0183/ 


c gk w is the planck-weighted k-distribution function due to h2o 

c line absorption in the 3 subbands (800-720,620-720,540-620 /cm) 

c of band 3 given in table 7. Note that the order of the sub-bands 

c is reversed. 


data gkw/ 0.1782,0.0593,0.0215,0.0068,0.0022,0.0000, 

2 0.0923,0.1675,0.0923,0.0187,0.0178,0.0000, 

3 0.0000,0.1083,0.1581,0.0455,0.0274,0.0041/ 


c ne is the number of terms used in each band to compute water vapor 

c continuum transmittance (table 6) . 

data ne /0,0,3, 1 , 1 , 1 ,0,0/ 


■tco2 are the six exp factors between levels kl and k2 

tran is the updated total transmittance between levels kl and k2 


c th2o is the 6 exp factors between levels kl and k2 due to 

c h2o line absorption. 


c tcon is the 3 exp factors between levels kl and k2 due to 

c h2o continuum absorption. 


■trnth2o is the total transmittance between levels kl and k2 due 
to both line and continuum absorption computed from eq. (52). 


do j-1 ,n 
do i=l,m 
th2o(i, j ,1) 
th2o(i, j ,2) 
th2o(i, j ,3) 
th2o(i, j ,4) 
th2o(i, j ,5) 
th2o(i, j ,6) 
enddo 
enddo 


th2o(i, j , 1) *h2oexp(i, j ,k, 1) 
th2o(i, j ,2)*h2oexp(i, j ,k,2) 
th2o(i, j ,3)*h2oexp(i,j ,k,3) 
th2o (i , j ,4) *h2oexp(i , j ,k,4) 
th2o(i, j ,5)*h2oexp(i, j ,k,5) 
th2o(i» j ,6)*h2oexp(i, j ,k,6) 


if (ne(ib) .eq.0) then 
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do j=l,n 
do i=i,m 


trnth2o =(fkw(i ,ib)*th2oCi, j , i) 

+ lkw(2,ib)*th2o(i, j ,2) 

+ fkw(3, ib)*th2o(i, j ,3) 

+ fkw(4»ib)*th2o(i, j ,4) 
+ fkw(5,ib)*th2o(i, j ,5) 

+ fkw(6, ib)*th2o(i, j ,6)) 

tranCi, j)=tranCi, j)+trnth2o 

enddo 

enddo 

elseif (ne(ib) . eq. 1) then 


do j = i,n 
do i=i,n 

tconCi, j ,1) 83 tconCi, j , l)*conexpCi, j ,k, 1) 

trnth2o =(f kw(i , ib)*th2o(i, j , 1) 

* + fkv(2 , ib)*th2o(i, j ,2) 

* + fkw(3,ib)*th2o(i, j ,3) 

‘ + fkwC4,ib)*th2o(i, j ,4) 

‘ + fkw(5,ib)«th2o(i, j ,6) 

1 + fkw(6,ib)*th2o(i, j ,6) )*tconCi, j , 1) 

tr an C i , j ) =tr an ( i , j ) *trnth2o 

enddo 

enddo 

else 

do j = l ,n 
do i=i,m 

tconCi, j , 1) = tconCi, j ,i)*conexpCi, j ,k, 1) 
tconCi, j ,2) = tconCi, j ,2)*conexpCi, j ,k,2) 
tconCi, j ,3)= tconCi, j ,3)*conexpCi, j ,k,3) 

trnth2o - C gkwCl , 1) *th2oCi, j , 1) 

+ gkwC2, l)*th2oCi, j,2) 

+ gkwC3, i)*th2o(i, j ,3) 

+ gkwC4,l)*th2oCi. j,4) 

+ gkwC5,l)*th2o(i, j,5) 

+ gkwC6, l)*th2oCi, j ,6) ) * tconCi, j,i) 
+ C gkw C 1 , 2) *th2o Ci , j , 1) 

+ gkwC2,2)*th2oCi, j,2) 
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+ gkw(3,2)*th2o(i, j.3) 

+ gk»(4,2)*th2o(i, j ,4) 

+ gkw(5 ,2) *th2o(i, j ,5) 

+ gkw(6,2)*th2o(i, j,6) ) * tcon(i,j,2) 

+ ( gkw (1,3) *th2o (i , j , 1) 

+ gkw(2,3)*th2o(i, j ,2) 

+ gkw(3,3)*th2o(i, j ,3) 

+ gkw(4,3)*th2o (i, j ,4) 

+ gk>(5, 3)*th2o(i, j ,5) 

+ gkw(6,3)*th2o(i,j,6) ) * tcon(i,j.3) 


tran(i, j)=tran(i, j)*trnth2o 


enddo 

enddo 

end if 


return 

end 


7.5 Subroutine C02KDIS 


c **** ************************************************ ****************** 
subroutine co2kdis (m,n,np ,k, co2exp,tco2 , tran) 

C * * * * * * # # # * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * 4c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

c compute co2 transmit tances between levels kl and k2 for min soundings 

c using the k-distribution method with linear pressure scaling, 
c 

c computations follow eq. (34) of Chou and Suarez (1994) . 
c 

c input parameters 

c number of grid intervals in zonal direction (m) 

c number of grid intervals in meridional direction (n) 

c 

c updated parameters 

c transmittance between levels kl and k2 due to co2 absorption 
c for the various values of the absorption coefficient (tco2) 
c total transmittance (tran) 
c 

c ************ ********************************************************** 
implicit none 
integer m,n,np,k,i,j 

c input parameters 

real*8 co2exp(m,n,np,6,2) 

c updated parameters 

real*8 tco2(m,n,6, 2) ,tran(m,n) 

c static data 

real*8 gkc(6,2) 

c temporary arrays 

real*8 xc 

c gkc is the planck-weighted co2 k-distribution function 

c in the band-wing and band-center regions given in table 7. 
c for computing efficiency, sub-bands 3a and 3c are combined. 

data gkc/ 0.1395,0.1407,0.1549,0.1357,0.0182,0,0220, 

2 0.0766,0.1372,0.1189,0.0335,0.0169,0.0059/ 

c tco2 is the 6 exp factors between levels kl and k2. 

c xc is the total co2 transmittance given by eq. (53) . 

do j=l,n 
do i=l,m 
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c 


band- wings 


tco2(i, j ,1 , l)=tco2(i, j , 1, i)*co2exp (i , j ,k,i,l) 
xc* gkc(l,l)*tco2(i, j ,1,1) 

tco2(i, j *2, l)=tco2(i f j ,2, i)*co2exp(i, j >k,2,l) 
xc=xc+gkc(2 , i)*tco2(i, j ,2,1) 

tco2(i, j , 3 , D=tco2(i, j,3,l)*co2exp(i, j ,k,3,l) 
xc=xc+gkc(3 , 1) *tco2(i , j ,3,1) 

tco2(i, j ,4, l)=tco2(i, j ,4, l)*co2exp(i, j ,k,4,l) 
xc=xc+gkc(4, l)*tco2(i, j ,4,1) 

tco2(i, j ,5 , l)*tco2(i, j ,5, 1) *co2exp(i , j ,k,5,l) 
xc=xc+gkc(5, l)*tco2(i, j ,5,1) 

tco2(i, j ,6 , l)=tco2(i, j ,6, l)*co2exp(i, j ,k,6,l) 
xc=xc+gkc(6, l)*tco2(i, j ,6,1) 

c band-center region 

tco2(i, j , 1 ,2)=tco2(i, j , 1 ,2)*co2exp(i, j ,k, 1 ,2) 
xc=xc+gkc(l ,2)*tco2(i, j ,1,2) 

tco2(i, j ,2,2)=tco2(i, j , 2 , 2) *co2exp(i , j ,k,2,2) 
xc=xc+gkc(2 ,2)*tco2(i, j ,2,2) 

tco2(i, j ,3 , 2)=tco2(i, j ,3,2)*co2exp(i, j ,k,3,2) 
xc=xc+gkc(3,2)*tco2(i, j ,3,2) 

tco2(i, j ,4,2)=tco2(i, j , 4,2)*co2exp(i, j ,k,4,2) 
xc=xc+gkc (4 ,2) *tco2(i , j ,4,2) 

tco2(i, j ,5 ,2)=tco2(i, j , 5 , 2) *co2exp(i , j ,k,5,2) 
xc=xc+gkc (5 ,2)*tco2(i , j ,5,2) 

tco2(i, j ,6,2)=tco2(i, j , 6,2) *co2exp(i , j ,k,6,2) 
xc=xc+gkc(6 ,2) *tco2(i, j ,6,2) 

t r an ( i , j ) ^tr an ( i , j ) *xc 

enddo 

enddo 

return 

end 
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7.6 Subroutine H20EXPS 


subroutine h2oerps (ib,m,n l np # dh2o ,pa,dt ,h2oexp) 
c**********^*********************************************************** 
c compute exponentials for water vapor line absorption 
c in individual layers . 
c 

c input parameters 

c spectral band (ib) 

c number of grid intervals in zonal direction (m) 
c number of grid intervals in meridional direction (n) 
c number of layers (np) 

c layer water vapor amount for line absorption (dh2o) 
c layer pressure (pa) 
c layer temperature minus 250K (dt) 
c 

c output parameters 

c 6 exponentials for each layer (h2oexp) 
c 

C**********************************************************************' 

implicit none 

integer ib,m,n,np, i, j ,k,ik 

c input parameters 

real*8 dh2o(m,n,np) ,pa(m,n,np) ,dt (m,n,np) 

c output parameters 

real*8 h2oexp(m,n,np, 6) 

c static data 

integer mw(8) 

real*8 xkw (8) , aw(8) ,bw(8) 

c temporary arrays 

real*8 xh 


c xkw are the absorption coefficients for the first 

c k-distribution function due to water vapor line absorption 
c (tables 4 and 7) . units are cm*+2/g 

data xkw / 29.55 , 4.167e-i, 1.328e-2, 5.250e-4, 

* 5 . 25e-4, 2.3406-3, 1.320e-0, 5.250e-4/ 

c are the ratios between neighboring absorption coefficients 

c for water vapor line absorption (tables 4 and 7) . 
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data mw /6, 6, 8, 6, 6, 8, 6, 16/ 


c aw bw (table 3) are the coefficients for temperature scaling 

c in eq. (25) . 

data aw/ 0.0021, 0.0140, 0.0167, 0.0302, 

* 0.0307, 0.0154, 0,0008, 0.0096/ 
data bw/ -l.Oie-5, 5.57e-5, 8.54e— 5, 2.96e-4, 

* 2 . 86e-4, 7.53e-5,-3.52e-6, 1.64e-5/ 


c note that the 3 sub-bands in band 3 use the same set of xkw, aw, 
c and bw. therefore, h2oexp for these sub-bands are identical. 
c ********** ************************************************************ 


do k=l ,np 
do j=l ,n 
do i=l ,m 

c z h i s the scaled water vapor amount for line absorption 

c computed from (27) . 

xh 52 dh2o(i, j ,k)*(pa(i, j ,k)*0 .002) 

1 * ( 1 . +(aw(ib)+bw(ib)* dt(i, j ,k))*dt (i, j ,k) ) 

c h2oexp is the water vapor transmittance of the layer (k2-l) 

c due to line absorption 

h2oexp(i , j ,k, 1) * exp(-xh*xkw(ib) ) 

enddo 

enddo 

enddo 

do ik=2,6 

if (mw(ib) . eq . 6) then 

do k=l,np 
do j=l,n 
do i=l,m 

ih = h2oexp(i, j ,k, ik-1) *h2oexp(i , j ,k, ik-1) 
h2oexp(i, j ,k, ik) = xh*xh*xh 
enddo 
enddo 
enddo 

elseif (mw(ib) . eq . 8) then 

do k=l ,np 
do j=l ,n 
do i=l,m 
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rh = h2oexp(i, j ,k, ik-i) *h2oexp(i , j ,k, ik-i) 
xh = xh*xh 

h2oexp(i f j f k,ik) = xh*xh 
enddo 
enddo 
enddo 

else 

do k=l,np 
do j=l,n 
do i=i,m 

xh = h2oexp(i, j,k,ik-l)*h2oexp(i, j ,k,ik-i) 
xh = xh*xh 
xh = xh*xh 

h2oexp(i, j ,k, ik) = xh*xh 
enddo 
enddo 
enddo 

end if 
enddo 


return 

end 
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7.7 Subroutine CONEXPS 


C** ******** ******************************* ******* ********* ************* 

subroutine c onexps(ib,m,n,np, dcont .conexp) 

c ******* ***************************************** ********************** 
c compute exponentials for continuum absorption in individual layers. 

c 

c input parameters 

C spectral band (ib) 

c number of grid intervals in zonal direction (m) 
c number of grid intervals in meridional direction (n) 
c number of layers (np) 

c layer scaled water vapor amount for continuum absorption (dcont) 
c 

c output parameters 

c 1 or 3 exponentials for each layer (conexp) 
c 

Q ************************************ ********************************** 

implicit none 

integer ib,m,n,np, i, j ,k,iq 

c input parameters 

real*8 dcont (m,n,np) 

c updated parameters 

real*8 conexp (m,n,np, 3) 

c static data 

integer ne(8) 
real*8 xke (8) 


xke are the absorption coefficients for the first 
k-distribution function due to water vapor continuum absorption 
(table 6). units are cm**2/g 

data xke / 0.00, 0.00, 27.40, 15.8, 

. 9.40, 7.75, 0.0, 0.0/ 


■ne is the number of terms in computing water vapor 
continuum transmittance (Table 6) . 
band 3 is divided into 3 sub-bands. 


data ne /0, 0,3, 1,1, 1,0,0/ 


do k=i,np 
do j=i,n 
do i»l,m 

conexp(i, j ,k, 1) = exp(-dcont(i, j,k)*xke(ib)) 
enddo 
enddo 
enddo 

if (ib . eq. 3) then 

c the absorption coefficients for sub-bands 3b (iq=2) and 3a (iq=3) 

c are, respectively, double and quadruple that for sub-band 3c (iq-1) 
c (table 6) . 

do iq=2,3 
do k=l,np 
do j*i,n 
do i=l,m 

conerpCi, j.k.iq) = conerp(i, j ,k, iq-1) *conexp(i, j ,k,iq-l) 
enddo 
enddo 
enddo 
enddo 

end if 

return 

end 
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7.8 Subroutine C02EXPS 


c * ******************************************* ************************** 
subroutine co2exps (m,n,np,dco2 , pa,dt ,co2exp) 


c 



c compute co2 exponentials for individual layers. 


c 

c input parameters 

c number of grid intervals in zonal direction (m) 

c number of grid intervals in meridional direction (n) 

c number of layers (np) 

c layer co2 amount (dco2) 

c layer pressure (pa) 
c layer temperature minus 250K (dt) 
c 

c output parameters 

c 6 exponentials for each layer (co2exp) 

c********************************************************************** 

implicit none 
integer m,n,np,i,j,k 


c input parameters 

real*8 dco2(m,n,np) ,pa(m,n,np) ,dt(m,n,np) 


c output parameters 

real*8 co2exp(m,n,np,6,2) 


c static data 

real*8 ikc(2) ,ac(2) ,bc(2) ,pn(2) ,prc(2) 


c temporary arrays 

real*8 xc 

c xkc is the absorption coefficients for the 

c first k-distribution function due to co2 (table 7). 
c units are l/(cm-atm) stp . 

data xkc/2 . 656e-5,2. 656e-3/ 

c parameters (table 3) for computing the scaled co2 amount 

c using (27) . 

data prc/ 300.0, 30.0/ 

data pm / 0.5, 0.85/ 

data ac / 0.0182, 0.0042/ 
data be /I . 07e-4 ,2 . 00e-5/ 
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c ************************************************************** ******** 


do k=l,np 
do j=l,n 
do i=l,m 

c compute the scaled co2 amount from eq. (27) for band-wings 

c (sub-bands 3a and 3c) . 

zc = dco2(i, j,k)*(pa(i, j,k)/prc(l))**pm(l) 

1 *(i.+(ac(l)+bc(l)*dt(i, j,k))*dt(i, j ,k)) 

c six exponential by powers of 8 (table 7) . 

co2exp(i, j ,k, 1, l)=exp(-xc*xkc (1) ) 

xc=co2exp(i, j ,k, i, l)*co2exp(i, j »k,l,l) 
xc-xc*xc 

co2exp(i, j ,k,2, l)=xc*xc 

xc=co2exp(i, j ,k,2, l)*co2exp(i, j ,k,2, i) 
xc=xc*xc 

co2exp(i, j ,k, 3, l)=xc*xc 

xc=co2exp(i, j ,k,3, l)*co2exp(i, j ,k,3, i) 
xc=xc*xc 

co2exp(i, j ,k,4, l)=xc*xc 

xc=co2exp(i, j ,k,4, i)*co2exp(i, j ,k,4, 1) 
xc=xc*xc 

co2exp(i, j ,k,5, l)=xc*xc 

xc=co2exp(i, j ,k,5, 1) *co2exp(i, j »k,5, 1) 
xc=xc*xc 

co2exp(i, j ,k,6, l)=xc*xc 

c compute the scaled co2 amount from eq. (27) for band-center 

c region (sub-band 3b) . 

xc = dco2(i, j ,k)*(pa(i, j ,k)/prc(2) )**pm(2) 

1 *(1 .+(ac(2)+bc(2)*dt (i, j ,k))*dt (i, j ,k)) 

co2exp(i, j , k, 1, 2)=exp(-xc*xkc(2) ) 

xc=co2exp(i, j ,k,l,2)*co2exp(i, j ,k, 1,2) 
xc=xc*xc 

co2exp(i, j ,k, 2, 2)=xc*xc 

xc=co2exp(i, j ,k,2,2)*co2exp(i, j ,k,2,2) 
xc=xc*xc 

co2exp(i, j ,k, 3, 2)=xc*xc 

xc=co2exp(i, j ,k,3,2)*co2exp(i, j ,k,3,2) 
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xc=xc*xc 

co2exp(i, j,k,4,2)*xc*xc 

xc=co2exp(i, j ,k,4,2)*co2exp(i, j ,k,4,2) 
xc=xc*xc 

co2exp(i, j ,k,5,2)=xc*xc 


xc=co2exp(i , j , k , 5 , 2) *co2exp(i , j ,k , B , 2) 

XC“XC*XC 

co2exp(i, j ,k,6,2)«xc*xc 

enddo 

enddo 

enddo 

return 

end 
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7.9 Sample Program 


c*************** source program lor CLIRAD IR: October 1994 ************ 

implicit none 

integer m,n,ndim,np, i, j ,k 

parameter (m=2,n=2,ndim=250,np=75) 

c input parameters 

real*8 pl(m,ndim,np+l) ,ta(m,ndim,np) ,wa(m,ndim,np) , 

* oa(m,ndim,np) , ts(m,ndim) ,taucl (m,ndim,np) , ccld(m,ndim,np) 
real+8 co2 


c output parameters 

real*8 f lx(m,ndim,np+l) ,f lc(m,ndim,np+i) ,dfdts(m,ndim,np+i) 
real*8 coolr (m,ndim,np) ,st4(m,ndim) 

real*8 rx 


logical high 
data high/. true./ 
c data high/ .false ./ 


c mid-latitude summer profiles 


data (pi (1,1 

, i) , i=l,np+l)/ .00, 

.0006244, 

.0008759, 

.001229, 

1 

.001723, 

.002417, 

.003391, 

.004757, 

.006672, 

.009359, 

2 

.01313, 

.01842, 

.02583, 

.03623, 

.05082, 

.07129, 

3 

0.10, 

0.14, 

0.20, 

0.28, 

0.39, 

0.54, 

4 

0.76, 

1.07, 

1.50, 

2.10, 

2.95, 

4.14, 

5 

5.80, 

8.14, 

11.42, 

16.01, 

22.46, 

31.51, 

6 

44.20, 

62.00, 

85.78, 

109.55, 

133.32, 

157.10, 

7 

180,88, 

204.65, 

228.43, 

252.20, 

275.98, 

299.75, 

8 

323.52, 

347.30, 

371,08, 

394.85, 

418.63, 

442.40, 

9 

466.17, 

489.95, 

513.72, 

537.50, 

561.28, 

585.05, 

a 

608.83, 

632.60, 

656.38, 

680.15, 

703.92, 

727.70, 

b 

751.47, 

775.25, 

799.03, 

822.80, 

846.58, 

870.35, 

c 

894.13, 

917.90, 

941.67, 

965.45, 

989.22, 

1013.00/ 

data (ta(l,l, 

, i) , i=l ,np)/209. 86, 

210.20, 

210.73, 


1 

211.27, 

211.81, 

212.35, 

212.89, 

213.44, 

213.98, 

2 

214.53, 

215.08, 

215.62, 

216.17, 

216.74, 

218.11, 

3 

223.20, 

230.04, 

237.14, 

244.46, 

252.00, 

259.76, 

4 

267.70, 

274.93, 

274.60, 

269.38, 

262.94, 

256.45, 

5 

250.12, 

244.31, 

238.96, 

233.74, 

228.69, 

224.59, 

6 

221.75, 

219.10, 

216.64, 

215.76, 

215.75, 

215.78, 

7 

216.22, 

219.15, 

223.79, 

228.29, 

232.45, 

236.33, 

8 

239.92, 

243.32, 

246.53, 

249.56, 

252.43, 

255.14, 

9 

257.69, 

260.11, 

262.39, 

264.57, 

266.66, 

268.67, 

a 

270.60, 

272.48, 

274.29, 

276.05, 

277.75, 

279.41, 
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b 281.02. 282.59, 284.09, 285.53, 286.86, 288.06, 

c 289.13, 290.11, 291.03, 291.91, 292.76, 293.59/ 

data (wa(l,l,i) , i=l,np)/ 16*0.400E-05. 

1 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 , 0 .400E-05 , 0 . 400E-05 , 

2 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-06 , 0 . 400E-05 , 0 . 400E-05 . 

3 0 . 400E-05 . 0 . 400E-05 , 0 . 400E-06 , 0 .400E-05 , 0 . 400E-05 ,0 . 400E-05 , 

4 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 , 0 . 400E-05 ,0 . 406E-06 , 

5 0 . 520E-05 ,0.1 15E-04 , 0 . 275E-04, 0 . 572E-04 , 0 . 107E-03 ,0 . 166E-03 . 

6 0 . 223E-03 , 0 . 285E-03 , 0 . 360E-03 , 0 . 446E-03 , 0 . 547E-03 ,0 . 655E-03 , 

7 0 . 767E-03 , 0 . 890E-03 , 0 . 103E-02 ,0.1 18E-02 , 0 . 136E-02 ,0 . 159E-02 , 

8 0 . 190E-02.0. 225E-02, 0. 264E-02.0 . 306E-02.0 . 351E-02.0 .399E-02 , 

9 0 . 450E-02 , 0 . 504E-02 , 0 . 560E-02 , 0 . 6 19E-02 , 0 . 680E-02 , 0 . 742E-02 , 

a 0 . 805E-02 , 0 . 869E-02 , 0 . 935E-02 , 0 . 100E-0 1,0. 107E-01 .0 . 1 13E-01/ 


c 


data (oa(l , 1 , i) , i=i ,np) / 


1 0 . 6427E-07, 

2 0 . 3779E-06 , 

3 0 . 6026E-06 , 

4 0 . 1009E-05, 

5 0 . 2716E-05 , 

6 0 . 7649E-05 , 

7 0 . 9900E-05 , 

8 0 . 2604E-0S , 

9 0 . 4483E-06, 
a 0. 1932E-06, 
b 0 . 1270E-06 , 
c 0 . 9238E-07 , 
d 0 . 7490E-07 , 
e 0 . 6258E-07 , 
f 0 . S517E-07 , 


0 . 2022E-06, 
0.4224E-06, 
0.6482E-06, 
0 . 1313E-05 , 
0.3U7E-05, 
0.9102E-05, 
0.8525E-05, 
0.1523E-05, 
0 . 3518E-06 , 
0. 1761E-06, 
0. 1180E-06, 
0 . 8828E-07 , 
0.7205E-07, 
0.6069E-07, 
0.5398E-07, 


0 . 2458E06 , 
0.4671E-06, 
0 . 6940E-06, 
0 . 1635E-05 , 
0 . 3590E-05 , 
0.9598E-05, 
0 . 7098E-05 , 
0. 1025E-05, 
0 . 3021E-06 , 
0 . 1S99E-06, 
0 . 1094E-06 , 
0 . 8458E-07 , 
0 . 6939E-07 , 
0 . 5928E-07 , 
0.5281E-07, 


0.2896E-06, 
0 . 5120E-06, 
0.740 IE-06, 
0 . 1976E-05, 
0.4645E-05, 
0 . 9944E-05 , 
0.5761E-05. 
0.7859E-06, 
0.2516E-06, 
0. 1466E-06, 
0. 1028E-06, 
0 . 8098E-07, 
0.6706E-07, 
0.5789E-07, 
0.5167E-07, 


0.3337E-06, 
0.5572E-06, 
0.7934E-06, 
0.2336E-05, 
0.5897E-05, 
0 . 1008E-04, 
0.4231E-05, 
0.5983E-06, 
0.2117E-06, 
0. 1366E-06, 
0.9748E-07, 
0.7782E-07, 
0.6480E-07, 
0.5651E-07, 
0.5054E-07/ 


ts(l . 1)=294 .0 

■m x n is the number of atmospheres 


do j = i,n 
do i=i,m 
ts(i , j)=ts(l , 1) 
enddo 
enddo 


■np is the number of atmospheric layers 


do k=i ,np 
do j=l,n 
do i=l ,m 

pl(i, j ,k)=pl(l,l,k) 
ta(i, j ,k)=ta(i , 1 ,k) 
sa(i, j ,k)=wa(l, l,k) 
oa(i, j ,k)=oa(l, i,k) 
enddo 
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enddo 

enddo 


np+i is the index for the surface level 

do j=i,n 
do i^l.m 

pl(i, j,np+i)=pl(l,l,np+l) 
enddo 
enddo 

assign co2 amount . units are parts/part by volumn 
co2=300 . e-6 

specify cloud cover and optical thickness 

do k=l,np 
do j=i,n 
do i=l,m 
ccld(i, j ,k)=0 .0 
taucl(i, j ,k)=0. 0 
enddo 
enddo 
enddo 

do k=46,49 
do j-i,n 
do i=i,m 
ccld(i, j ,k)=0. 5 
taucl (i, j ,k)=2 . 5 
enddo 
enddo 
enddo 

compute ir fluxes due to vater vapor, co2, and o3 

call irrad (m,n, ndim,np, taucl, ccld, pi, ta,wa, 

oa,co2 ,ts ,high,flx,f lc,dfdts,st4) 

compute cooling rate profile 

do k=l,np 
do j=i,n 
do i=i,m 

coolrCi, j ,k)=(flx(i, j,k+l)-f lx (i, j,k) >*8.441874/ 
(pl(i, j ,k+l)-pl(i, j,k)) 

enddo 

enddo 

enddo 

print out input data 


writ® (7,100) 

writ® (7,102) ts(l , 1) ,co2 

writ® (7,10B) 

writ® (7,108) 

writ® (7,110) 


do j-1,1 
do i-1,1 
do k=l,np 

writ® (7,115) k,pl(i,j,k), tlx (i, j ,k) ,flc(i, j ,k) ,dfdts(i, j ,k) 
writ® (7,116) ta(i, j ,k) ,wa(i, j ,k) ,oa(i, j ,k) ,ccld(i, j ,k) , 

* taucl(i, j ,k) , coolr(i,j,k) 


enddo 

enddo 

enddo 


i-1 

j-i 

k=np+l 

write (7,116) k,pl(i,j.k), f lx(i, j ,k) ,flc(i, j ,k) ,dfdts(i,j ,k) 
writ® (7,117) st4(i,j) 


100 format (/,* **♦* Sample mid-latitude summer atmosphere »***», /) 
102 format (/,' ts= • , f 8 . 2, ' K > ,/ , ’ co2= ’ ,e9. 3, ’ ppmw’ ./) 

105 format (/,' ****** INPUT DATA ***♦**>, 28x, 

* ******* RESULTS ******»,/) 

108 format (lOx, 'p’ ,6x, *T* ,7*. ’q’ ,8x, >o3> ,6r, ’cold’ ,2x, ’taucl’ . 

* 3x, ’fix’ ,6x, ’flc’ ,4x, ’dfdts’ ,2x, ’cooling rate’) 

110 format (8x, ’ (mb) * ,4x, ’ (k) ’ ,4x, ’ (g/g) 1 ,6x, ’ (g/g) ’ ,16x, ’ (W/m‘2) ’ . 

* lx, ’ (V/m*2) ’ lx, ’ (W/m"2/K) ' , lx, 1 (C/day) ’ ,/) 

115 format (i4,fl0.4, 40x,f 7 . 2,2x,f7 . 2, lx,f6 .2) 

116 format (14x,f7 . 2,2el0 .3,f 5 . 2,f 6 .2, 25xf8.2) 

117 format (/,2x, ’st4=’ ,f8.2, ’ W/m**2’) 

stop 

end 
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7.10 Verification Output of Sample Program 


**** Sample mid-latitude summer atmosphere **** 


ts= 294.00 K 
co2=0.300E-03 ppmv 


IIPUT DATA ****** ****** RESULTS ****** 

P T q o3 ccld taucl fix flc dfdts cooling rate 

(mb) (k) (g/g) (g/g) (W/m"2) (W/m-2) (W/m'2/I) (C/day) 


1 

0.0000 



-191.22 

-293.07 

-0.08 




209,86 0.400E-05 0.643E-07 0.00 

0.00 




91.63 

2 

0 . 0006 



-191.21 

-293.07 

-0.08 




210.20 0.400E-05 0.202E-06 0.00 

0.00 




24.69 

3 

0.0009 



-191.21 

-293.07 

-0.08 




210.73 0.400E-05 0.246E-06 0.00 

0.00 




32.02 

4 

0.0012 



-191.21 

-293.06 

-0.08 




211.27 0. 400E-05 0.290E-06 0.00 

0.00 




22.13 

5 

0.0017 



-191.21 

-293.06 

-0.08 




211.81 0.400E-05 0.334E-06 0.00 

0.00 




22.57 

6 

0.0024 



-191.20 

-293.06 

-0.08 




212.35 0.400E-05 0.378E-06 0.00 

0.00 




15.31 

7 

0.0034 



-191.20 

-293.06 

-0.08 




212.89 0 . 400E-05 0.422E-06 0.00 

0.00 




14.91 

8 

0.0048 



-191.20 

-293.06 

-0.08 




213.44 0 . 400E-05 0.467E-06 0.00 

0.00 




9.93 

9 

0.0067 



-191.20 

-293.05 

-0.08 




213.98 0. 400E-05 0.512E-06 0.00 

0.00 




7.08 

10 

0.0094 



-191.19 

-293.05 

-0.08 




214.53 0.400E-05 0.557E-06 0.00 

0.00 




7.29 

11 

0.0131 



-191.19 

-293.05 

-0.08 




215.08 0. 400E-05 0.603E-06 0.00 

0.00 




5.53 

12 

0.0184 



-191.19 

-293.05 

-0.08 




215.62 0.400E-05 0.648E-06 0.00 

0.00 




4.28 

13 

0.0258 



-191.18 

-293.04 

-0.08 




216.17 0.400E-05 0.694E-06 0.00 

0.00 




3.13 

14 

0.0362 



-191.18 

-293.04 

-0.08 




216.74 0.400E-05 0.740E-06 0.00 

0.00 




1.93 

15 

0.0508 



-191.18 

-293.03 

-0.08 




218.11 0.400E-05 0.793E-06 0.00 

0.00 




0.28 

16 

0.0713 



-191.18 

-293.03 

-0.08 




223.20 0 . 400E-05 0.101E-05 0.00 

0.00 




1.50 

17 

0.1000 



-191.17 

-293.03 

-0.08 




230.04 0. 400E-05 0.131E-05 0.00 

0.00 




2.19 

18 

0.1400 



-191.16 

-293.02 

-0.08 




237.14 0.400E-05 0.164E-05 0.00 

0.00 




3.12 

19 

0.2000 



-191.14 

-293.00 

-0.08 




244.46 0. 400E-05 0.198E-05 0.00 

0.00 




3.88 

20 

0.2800 



-191.10 

-292 .96 

-0.08 




252.00 0.400E-05 0.234E-05 0.00 

0.00 




5.24 

21 

0.3900 



-191.03 

-292.89 

-0.08 




259.76 0 . 400E-05 0.272E-05 0.00 

0.00 




7.18 

22 

0.5400 



-190.91 

-292.77 

-0.08 

' ' ' 



267.70 0.400E-05 0.312E-05 0.00 

0.00 




8.92 
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23 

0.7600 


-190.67 

-292.54 -0.08 



274.93 0. 400E-05 0.359E-05 0.00 

0.00 



12.15 

24 

1,0700 


-190.23 

-292.10 -0.08 



274.60 0.400E-05 0.465E-05 0.00 

0.00 



11.64 

25 

1.5000 


-189.63 

-291.51 -0.08 



269.38 0.400E-05 0.590E-05 0.00 

0.00 



9.29 

26 

2.1000 


-188.97 

-290.87 -0.08 



262.94 0.400E-05 0.765E-05 0.00 

0.00 



7.26 

27 

2.9500 


-188.24 

-290.16 -0.08 



256.45 0 .400E-05 0.910E-05 0.00 

0.00 



5.72 

26 

4.1400 


-187.44 

-289.39 -0.08 



250.12 0.400E-05 0.960E-05 0.00 

0,00 



4.85 

29 

5.8000 


-186.48 

-288.48 -0.08 



244.31 0. 400E-05 0.994E-05 0.00 

0.00 



4.03 

30 

8,1400 


-185.37 

-287,44 -0.08 



238.96 0. 400E-05 0.101E-04 0.00 

0.00 



3.19 

31 

11.4200 


-184.13 

-286.30 -0.08 



233.74 0.400E-05 0.990E-05 0.00 

0.00 



2.56 

32 

16.0100 


-182.73 

-285.04 -0.08 



228.69 0.400E-05 0.853E-05 0.00 

0.00 



2.24 

33 

22.4600 


-181.02 

-283.52 -0.08 

1.51 


224.59 0.400E-05 0.710E-05 0.00 

0.00 



34 

31.5100 


-179,40 

-282.17 -0.08 

1.27 


221.75 0.400E-05 0.576E-0S 0.00 

0,00 



35 

44.2000 


-177.50 

-280,66 -0.08 

1.00 


219.10 0 . 400E-05 0.423E-05 0.00 

0.00 



36 

62.0000 


-175.40 

-279.08 -0.08 

0.58 


216.64 0.400E-05 0.260E-05 0.00 

0.00 



37 

85.7800 


-173.75 

-278.09 -0.08 

0.43 


215.76 0.400E-05 0.152E-05 0.00 

0.00 



38 

109.5500 


-172.54 

-277.40 -0,08 

0.39 


215.75 0. 400E-05 0.102E-05 0.00 

0.00 



39 

133.3200 


-171.45 

-276.75 -0.08 

0.28 


215.78 0.406E-05 0.786E-06 0.00 

0.00 



40 

157.1000 


-170.66 

-276.39 -0.08 

0.10 


216.22 0.520E-05 0.598E-06 0.00 

0.00 



41 

180.8800 


-170.38 

-276.51 -0.09 

0,48 


219.15 0.115E-04 0.448E-06 0.00 

0.00 



42 

204.6500 


-169.02 

-275.53 -0.09 

1.16 


223.79 0. 275E-04 0.352E-06 0.00 

0.00 



43 

228.4300 


-165.76 

-272.66 -0.09 

1.80 


228.29 0 . 572E-04 0.302E-06 0,00 

0.00 



44 

252.2000 


-160.69 

-268.06 -0.09 

2.18 


232.45 0.107E-03 0.252E-06 0.00 

0.00 



45 

275.9800 


-154.55 

-262.55 -0.09 

2.32 


236.33 0.166E-03 0.212E-06 0.00 

0.00 



46 

299.7500 


-148.02 

-256.92 -0.09 

19.25 


239.92 0.223E-03 0.193E-06 0.50 

2.50 



47 

323.5200 


-93.83 

-251.64 -0.17 

6.42 


243.32 0.285E-03 0.176E-06 0.50 

2.50 



48 

347.3000 


-75.74 

-246.54 -0.34 

-2.61 


246.53 0.360E-03 0.160E-06 0.50 

2.50 



49 

371.0800 


-83.10 

-241.45 -0.67 

-12.85 


249.56 0. 446E-03 0.147E-06 0.50 

2.50 



50 

394.8500 


-119.28 

-236.35 -1.31 

-0.15 


252.43 0 . 547E-03 0.137E-06 0.00 

0.00 



51 

418.6300 


-119.71 

-231.24 -1.32 

0.02 


255.14 0.655E-03 0.127E-06 0.00 

0.00 



52 

442.4000 


-119.66 

-226.19 -1.32 

0.12 


257.69 0.767E-03 0.118E-06 0.00 

0.00 



53 

466.1700 


-119.31 

-221.29 -1.32 
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260.11 0 . 890E-03 0.109E-06 0.00 

0.00 



0.22 

64 

489.9500 


-118.68 

-216.44 -1.33 


262.39 0. 103E-02 0.103E-06 0.00 

0.00 



0.29 

55 

513.7200 


-117.88 

-211.61 -1.33 


264,57 0.118E-02 0.975E-07 0.00 

0.00 



0.36 

56 

537 . 5000 


-116.87 

-206.82 -1.33 


266.66 0. 136E-02 0.924E-07 0.00 

0.00 



0.42 

67 

561.2800 


-115.67 

-201.99 -1.34 


268.67 0 . 159E-02 0.883E-07 0.00 

0.00 



0.50 

58 

585.0500 


-114.26 

-196.97 -1.34 


270.60 0 . 190E-02 0.846E-07 0.00 

0.00 



0.56 

59 

608.8300 


-112.70 

-191.82 -1.35 



272.48 0.225E-02 0.810E-07 0.00 

0.00 



0.63 

60 

632.6000 


-110,93 

-186.36 -1.36 



274.29 0 . 264E-02 0.778E-07 0.00 

0.00 



0.67 

61 

656.3800 


-109.03 

-180.95 -1.37 



276.05 0.306E-02 0.749E-07 0.00 

0.00 



0.74 

62 

680.1500 


-106.95 

-175.29 -1.38 


277.75 0.351E-02 0.721E-07 0.00 

0.00 



0.79 

63 

703.9200 


-104.72 

-169.63 -1.39 


279.41 0.399E-02 0.694E-07 0.00 

0.00 



0.86 

64 

727.7000 


-102.29 

-163.82 -1.41 


281.02 0.450E-02 0.671E-07 0.00 

0.00 



0.92 

65 

751.4700 


-99.69 

-157.89 -1.44 


282.59 0. 504E-02 0.648E-07 0.00 

0.00 



1.02 

66 

775.2500 


-96.81 

-151.75 -1.46 



284.09 0.560E-02 0.626E-07 0.00 

0.00 



1.09 

67 

799.0300 


-93.74 

-145.49 -1.50 



285.53 0.619E-02 0.607E-07 0.00 

0.00 



1.21 

68 

822.8000 


-90.34 

-139.01 -1.54 



286.86 0.680E-02 0.593E-07 0.00 

0.00 



1.31 

69 

846.5800 


-86.64 

-132.32 -1.59 



288.06 0 . 742E-02 0.579E-07 0.00 

0.00 



1.38 

70 

870.3500 


-82.75 

-125.49 -1.66 



289.13 0.805E-02 0.565E-07 0.00 

0.00 



1.41 

71 

894.1300 


-78.77 

-118.64 -1.75 



290.11 0. 869E-02 0.552E-07 0.00 

0.00 



1.44 

72 

917.9000 


-74.72 

-111.85 -1.86 


291.03 0.935E-02 0.540E-07 0.00 

0.00 



1.50 

73 

941.6700 


-70.49 

-104.80 -2.01 



291.91 0 . 100E-01 0. 528E-07 0.00 

0.00 



1.58 

74 

965.4500 


-66.03 

-97.57 -2.24 



292.76 0. 107E-01 0.517E-07 0.00 

0.00 



1.70 

75 

989.2200 


-61.23 

-90.04 -2.64 



293.59 0.113E-01 0.505E-07 0.00 

0.00 



2.22 

76 1013.0000 


-54.99 

-81.15 -5.76 



st 4= -423.62 W/m*»2 
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